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Orthogonal  wavelet  transforms  have  been  applied  to  the  field  of  signal  and 
image  processing  with  promising  results  in  compression  and  denoising.  Coefficients  of 
such  a  transform  constitute  a  complete  representation  of  a  signal  without  redundancy. 
However,  there  are  applications  where  complete  representations  are  disadvantageous. 
In  this  thesis,  we  examine  classes  of  Fourier  transforms  and  wavelet  transforms  in 
terms  of  their  efficacy  of  representing  convolution  operators.  We  have  identified  two 
shortcomings  associated  with  complete  representations  of  the  discrete-time  domain: 
(1)  the  lack  of  translation  invariance  and  the  (2)  a  possible  anomaly  of  aliasing- 
enhancement. 

On  the  other  hand,  our  analysis  showed  that  overcomplete  wavelet  represen- 
tations do  not  bear  those  shortcomings  of  their  non-redundant  counterparts.  Our 
framework  of  overcomplete  wavelet  representations  include  construction  algorithms 
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and  prototype  filters,  spatial-frequency  interpretation  and  three  operations.  Capabil- 
ities of  spatial-frequency  localization  were  quantitatively  evaluated  using  uncertainty 
factors.  Associated  with  gain,  shrinking  and  envelope  operators,  algorithms  for  con- 
volution, denoising  and  analysis  of  power  density  distribution  were  presented  and 
analyzed. 

The  framework  of  overcomplete  wavelet  representations  was  then  applied  to 
segmentation  of  textured  images  and  image  deblurring.  We  demonstrated  that  en- 
velopes as  feature  vectors  performed  well  in  segmenting  both  natural  and  synthetic 
textures.  We  showed  that  gain  and  shrinking  operators  may  be  used  for  image  de- 
blurring  and  discuss  limitations  of  the  mothodology. 
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CHAPTER  1 
INTRODUCTION 

1.1    Motivations  of  the  Research 

In  the  past  decade,  wavelets  have  stimulated  great  interests  in  many  fields 
of  applied  mathematics  and  engineering,  generating  a  significant  amount  of  research 
activities.  Within  these  fields,  traditional  theories  and  techniques  have  been  reeval- 
uated to  explore  opportunities  for  wavelet  applications. 

Although  continuous  time  wavelets  possess  mathematical  elegance,  discrete- 
time  wavelet  transforms  are  of  special  importance  for  practical  applications.  Among 
them,  discrete-time  orthogonal  wavelet  transforms  have  gained  popularity.  Coeffi- 
cients of  such  a  transform  constitute  a  complete  representation  of  a  signal  without 
redundancy.  Complete  representations  have  produced  promising  results  in  image 
compression  and  denoising.  However,  there  are  applications  where  complete  repre- 
sentations are  disadvantageous.  Understanding  the  strength  and  weakness  of  various 
representations  and  the  relationships  among  them  is  essential  for  satisfying  intel- 
lectual inquires  and  developing  new  methodologies.  In  this  thesis,  we  considered 
both  classes  of  Fourier-based  representations  and  wavelet-based  representations  with 
emphasis  in  the  discrete-time  domain.  We  found  that  the  overcomplete  wavelet  rep- 
resentations served  to  bridge  the  two  classes. 

The  merit  of  a  particular  representation  should  be  judged  by  its  capability  to 
facilitate  solving  problems.  Two  important  but  difficult  problems,  segmentation  of 
textured  images  and  image  deblurring,  were  selected  as  benchmarks.  Accordingly, 
we  value  a  representation  by  its  capability  to  extract  certain  features  from  signals, 
efficacy  of  representing  an  operation  and  invariance  under  certain  operations. 
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With  the  understanding  of  various  representations  and  the  problems  in  hand, 
we  developed  a  framework  of  discrete-time  overcomplete  wavelet  representation.  The 
framework  included  construction  algorithms  and  prototype  filters,  time-frequency  in- 
terpretation and  three  operations.  Algorithms  associated  with  gain  operators,  shrink- 
ing operators  and  envelope  detectors  were  developed.  The  framework  of  overcomplete 
wavelet  representations  was  then  applied  to  segmentation  of  textured  images  and  im- 
age deblurring.  This  demonstrated  that  envelopes  as  feature  vectors  performed  well 
in  segmenting  both  natural  and  synthetic  textures.  For  the  deblurring  application, 
our  formulation  extended  the  wavelet- vaguelette  inversion  into  discrete-time  domain 
and  was  applicable  to  more  blurring  sources. 

1.2    Review  of  Related  Work 

On  the  general  topic  of  wavelet  theory,  a  brief  history  and  comprehensive  re- 
view were  included  in  Jawerth  and  Sweldens'  survey  [32].  The  papers  by  Daubechies 
[11]  and  Mallat  [42,  41]  were  considered  influential  for  the  recent  developments.  The 
paper  by  another  wavelet  pioneer  Strang  [59]  was  also  explanatory.  The  orthonormal 
wave  packet  bases  is  a  generalization  of  orthogonal  wavelets,  was  introduced  by  Coif- 
man  and  Meyer  [10].  The  book  by  Daubechies  [13]  included  in-depth  mathematical 
discussions  and  real  design  examples.  The  book  by  Vaidyanathan  [63],  which  was 
written  for  signal  processing  audience,  provided  much  insight  into  the  relationship 
between  wavelets  and  multirate  filter  banks. 

The  lack  of  translation  invariance  of  discrete-time  orthogonal  wavelets  was  rec- 
ognized as  a  major  shortcoming  for  pattern  recognition  by  many  researchers  [41,  59]. 
Various  proposals  were  made  to  address  the  issue.  Multiscale  edge  representations 
were  investigated  by  Mallat  and  Zhong  [44],  Mallat  and  Hwang  [43],  and  Berman  and 
Baras  [5].  Steerable  filters  with  shiftabilities  were  studied  by  Simoncelli  et  al.  [58]. 
A  orthogonal  representation  with  invariance  of  the  subband  structure  was  proposed 
by  Pesquet  et  al.  [52]. 
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Using  energy  distribution  in  subbands  of  wavelet  packet  decomposition  as 
signatures  for  texture  classification  was  presented  by  Laine  and  Fan  [36].  A  similar 
approach  was  also  reported  by  Chang  and  Kuo  [8].  Earlier  results  using  envelopes  of 
overcomplete  wavelet  packet  representation  for  texture  segmentation  was  reported  by 
Laine  and  Fan  [35].  A  simplified  version  of  Chapter  5  of  this  thesis  was  later  published 
[37].  The  overcomplete  wavelet  representation  was  also  used  by  User  [62]  with  an  ad 
hoc  feature  extraction  algorithm.  A  methodology  combining  Markov  random  field 
and  orthogonal  wavelet  packet  transform  was  described  by  Bello  [4]. 

The  nonlinear  wavelet  shrinkage  for  denoising  was  due  to  Donoho  and  John- 
stone [17].  The  approach  was  applied  to  power  spectrum  estimation  [48],  and  to 
medical  image  processing  [1].  Nonlinear  denoising  based  on  wavelet  maxima  rep- 
resentation by  Mallat  and  Hwang  shared  some  similar  ideas  [43],  although  tracing 
singularities  through  the  scale  space  turned  out  to  be  much  more  difficult,  especially 
in  a  two-dimensional  space. 

Manipulating  wavelet  representations  to  enhance  certain  features  is  obviously 
a  good  idea.  Jawerth  et  al.  [31]  reported  image  enhancement  using  weight  factors  to 
modify  coefficients  of  orthogonal  wavelet  transforms.  Laine  et  al.  [38]  demonstrated 
various  approaches  for  enhancement  of  digital  mammographic  images.  Fan  and  Laine 
[20]  pointed  out  the  connection  between  linear  gain  enhancement  using  overcomplete 
dyadic  wavelet  representations  and  the  unsharp  masking,  and  introduced  a  nonlinear 
function  for  contrast  enhancement.  However,  due  to  the  lack  of  an  objective  definition 
of  image  quality,  those  enhancement  algorithms  were  mostly  ad  hoc  and  hard  to 
evaluate.  On  the  other  hand,  image  deblurring  is  considered  to  be  a  well  defined 
problem.  The  approach  adapted  by  Banham  et  al.  [3]  is  to  decompose  the  linear 
minimum  mean  square  error  (LMMSE)  filters  into  the  orthogonal  transform  domain 
with  advantage  of  reduced  alliance  on  the  assumption  of  global  stationarity.  However, 
their  proposal  was  basically  the  linear  filtering  and  no  optimal  subband  selection  was 
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considered.  The  mathematical  framework  of  the  wavelet- vaguelette  inversion  was  due 
to  Donoho  [14,  16].  The  limitations  were  that  it  was  derived  on  the  continuous  time 
domain  and  inapplicable  to  some  common  sources  such  as  Gaussian  and  uniform 
blurs. 

1.3    Overview  of  the  Thesis 

This  thesis  is  organized  as  follows: 

Chapter  2  reviews  signal  representations  derived  from  the  Fourier  transform 
and  the  wavelet  transform  in  the  context  of  being  operated  upon  by  a  convolution 
operator.  The  scrutiny  reveals  that  complete  representations  suffer  from  not  being 
translation  invariant  and  the  existence  of  the  aliasing-enhancement  anomaly,  and  are 
thus  inappropriate  for  certain  applications. 

Chapter  3  presents  the  framework  of  the  overcomplete  wavelet  representation 
with  two  particular  building  algorithm  and  filters.  The  time-frequency  localization 
property  of  overcomplete  wavelet  representations  is  examined  under  the  uncertainty 
principle. 

Chapter  4  introduces  gain,  shrinking  and  envelope  operators  applicable  upon 
an  overcomplete  wavelet  representation.  Algorithms  for  determining  specific  param- 
eters are  detailed  with  examples  and  comparisons  shown.  Both  Chapter  5  and  6 
are  applications  of  the  theory  of  overcomplete  wavelet  representations  developed  in 
previous  chapters. 

Chapter  5  deals  with  the  topic  of  image  texture  segmentation.  The  segmen- 
tation problem  is  formulated  in  the  paradigm  of  spatial-frequency  analysis,  and  a 
texture  feature  based  on  the  concept  of  power  density  distribution  is  proposed. 

Chapter  6  addresses  the  problem  of  image  deblurring  using  the  gain  operator 
for  deconvolution  and  the  shrinking  operator  for  nonlinear  denoising. 

Chapter  7  concludes  the  research  and  proposes  future  directions. 


CHAPTER  2 

SIGNAL  REPRESENTATIONS  AND  CONVOLUTION  OPERATORS 

2.1  Introduction 

For  a  given  signal  x,  it  is  often  beneficial  to  transform  (map)  x  into  another 
form  (domain),  denoted  x  X.  X  is  called  the  representation  of  x  in  the  transform 
domain.  Instead  of  directly  dealing  with  the  original  signal  x,  we  may  deal  with 
its  representation  X.  By  doing  so  we  usually  transform  a  particular  problem  into 
a  different  form  and  make  it  easier  to  solve.  Clearly,  different  applications  demand 
different  representations. 

The  applications  considered  in  this  thesis  are  texture  image  segmentation  and 
image  deblurring.  Accordingly,  the  time-invariability  and  the  efficacy  of  represen- 
tations are  the  most  important  factors.  Both  properties  may  be  revealed  by  their 
behavior  under  a  convolution  operator.  With  these  in  mind,  we  shall  review  Fourier 
transforms  and  wavelet  transforms  of  one  dimensional  signals.  To  facilitate  our  dis- 
cussions, we  used  the  following  basic  definitions: 

Definition  2.1.1  (Translation  Operator)  For  a  continuous  function  x(t),  a  trans- 
lation operator  is  defined  as  Ts  [x(t)]  —  x(t  —  s);  For  a  discrete  function  x(n),  a  trans- 
lation operator  is  defined  as  Td  [x(n)]  =  x{n  —  d). 

Definition  2.1.2  (Convolution  Operator)  For  a  continuous  function  f(t),  a  con- 
volution operator  is  defined  as  Tf  [x(t)]  =  x(t)  *  f(t)  =  P^:c(t)/(£  -  r)rfr;  For  a 
discrete  function  x(n),  a  convolution  operator  is  defined  as  Tf  [x(n)]  =  x(n)  *  f(n)  = 
T,m=-oox(m)f(n~rn)-  The  functions  f  (t)  and  f{n)  are  called  the  kernel.  Notice  that 
a  translation  operator  is  a  special  convolution  operator  with  kernel  f(t)  =  5(t  -  r) 
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(f(n)  =  S(n  -  d)).  Therefore,  we  will  use  T  for  both  translation  and  convolution 
operators. 

Definition  2.1.3  (Eigenfunction  and  Eigenvalue)  For  a  linear  operator  T ,  if 
there  exists  a  function  e,  such  that  T  [e]  =  Xe,  where  A  is  a  constant,  the  function  e 
is  called  the  eigenfunction  of  T,  and  the  value  A  is  called  the  eigenvalue. 

Definition  2.1.4  (Translation  invariant)  For  a  given  signal  x,  a  translation  op- 
erator T  and  a  mapping  operator  M.,  if  the  two  operators  are  commutable  such  that 

M[T(x)]  =  T[M(x)], 

the  mapping  is  called  translation  invariant,  and  the  representation  X  is  said  to  be  a 
translation  invariant  representation. 

2.2    Fourier  Transforms 

Frequency-domain  representation  by  Fourier  transform  played  a  major  role 
in  signal  analysis  and  linear  system  theory  [51,  49,  63].  The  efficacy  of  Fourier 
representation  is  due  to  the  fact  that  its  basis  functions  of  complex  exponentials 
are  eigenfunctions  of  a  convolution  operator  [49,  page  39].  For  a  linear  convolution 
operator  T,  the  eigenequation  is: 

Tf  [e*"]  =  F(w)e*", 

where  eigenvalue  F(u)  is  the  Fourier  transform  of  f(t).  An  almost  identical  form 
can  be  derived  for  the  discrete-time  Fourier  transform.  Fourier  representations  are 
particularly  powerful  for  signal  deblurring  since  it  transforms  a  deconvoltuion  problem 
into  the  simple  form  of  divisions. 

However,  Fourier  representations  may  not  be  good  at  capturing  frequency 
evolution  over  time.  As  time  variable  is  added  to  Fourier  transforms,  basis  functions  of 
short-time  Fourier  transforms  are  no  longer  eigenfunctions  of  a  convolution  operator. 


2.2.1  Fourier  Transform  (FT) 

The  Fourier  transform  pair  may  be  written  as: 

/oo 
x(t)e~]u,tdt  (FT), 
-oo 
/•oo 

x(t)  =  —  /    X(u)e]u3tdu   (Inverse  FT). 
A  linear  convolution  operator  Tf  can  be  characterized  by 

/oo  r         ,  roo 

X(u)Tf  ejuJt  doo=  X(u)F(u))ejultduj, 
-00  J —  oc 

or, 

Y(u>)  =  X{w)F(u>). 

This  shows  that  Fourier  representation  of  the  convolution  operator  T  is  the 
Fourier  transform  of  its  kernel.  The  beauty  of  this  representation  is  that  on  the 
Fourier  domain,  a  convolution  (deconvolution)  is  simply  multiplications  (divisions) 
of  the  both  representations. 

For  a  translation  operator  with  kernel  f(t)  =  5(t  —  s),  F(u)  =  e~iws ,  and  there- 
fore |F(w)|  =  |Ar(o;)|,  which  means  that  magnitude  of  Fourier  transform  is  translation 
invariant. 

2.2.2  Discrete-time  Fourier  Transform  (DTFT) 

The  discrete-time  Fourier  transform  may  be  written  as  [49]: 

oo 

X(e?w)  =  Yl  x(n)e-]u,n  (DTFT), 

n=  —  oo 

x(n)  =  —  P  X(jw)e?Mdw   (Inverse  DTFT). 
Similar  to  its  continuous  counterpart, 

y(n)  =  Tf  [x(n)}  =  i-  f  X{^)F{e?u)^du), 

or, 

Y(eju})  =  X(ejuJ)F(eJul). 

For  a  translation  operator  with  f(n)  =  6(n  -  d).  F(ejuJ)  =  e_jW  and  thus 
\Y(e^)\  =  \X(e^)\. 
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2.2.3    Short-time  Fourier  Transform  (STFT) 

In  order  to  capture  frequency  evolution  of  a  non-stationary  signal,  a  window 
w(t)  is  introduced  into  the  Fourier  transform.  The  short-time  Fourier  transform  may 
be  written  as: 

/OO 
x(r)w(t  -  T)e~]UJTdT  (STFT). 
-oc 

roc 

x(t)  =  - — —  /     X(u,  t)eJU'tdu)      (Inverse  STFT). 
2irw(0)  J-oc 

where  we  assumed  w(0)  ^  0. 

For  a  convolution  operator  Tj, 

Tf  [w(t  -  r)e^']  =  F(u,  t  -  r)ejut, 

where  F(u;,t)  is  the  STFT  representation  of  the  kernel  f(t).  It  means  that  in  general 
the  windowed  complex  exponentials  are  no  longer  eigenfunctions  of  a  convolution  op- 
erator. Consequently,  STFT  representation  of  y(t)  =  Tf  [x(t)]  is  another  convolution 
of  either 


or. 


/oo 
F(uJ  -T)x(r)e-^TdT, 
-oo 

/oo 
X(u,t-T)f(r)e-^TdT. 
-oo 


For  a  translation  operator  with  f(t)  =  S(t  —  s),  Y(u,t)  =  X(u,t  —  s)e~^s, 
and  thus  |F(u;.t)|  =  \X(u,t  —  s)\.  Therefore,  magnitude  of  STFT  transform  is  a 
time-invariant  representation. 

2.2.4    Discrete-time  Short-time  Fourier  Transform  (DTSTFT) 

The  discrete-time  short-time  Fourier  transform  may  be  written  as  [49,  53] 

oo 

X(ejw,n)=  J2  x(k)w(n-  k)e-Jujk,   (DTSTFT)  (2.1) 

k=— oo 

x(n)  =  — F  X(ejuJ,n)e^ndu).  (Inverse  DTSTFT)  (2.2) 
2irw(Q)  J-TT 


For  a  discrete  convolution  operator  T)  and  y(n)  =  Tj  [x(n)],  we  have 

oo 

Y{e?u,n)=  £  X(e^,n-k)f(k)e-^k, 

k=  —  oo 

or, 

00 

Y(ejuJ,n)  =  Y,  F{e]UJ,n-k)x{k)e-ju}k. 

k=  — oo 

For  a  translation  operator  7)  with  f(n)  =  8(n-k),  \Y(ejuJ,n)\  =  |A"(eJa,,n  -  /c)|. 
Thus,  the  magnitude  of  DTSTFT  is  also  a  translation-invariant  representation. 

2.2.5    Connection  Between  Short-time  Fourier  Transform  and  Filter  Banks 

The  discrete-time  short-time  Fourier  transform  X(e^,  n)  is  a  two-dimensional 
function  of  a  continuous  variable  uj  and  a  integer  variable  n.  Is  it  possible  to  have 
only  finite  samples  X{e?Uk,n)  in  the  frequency  domain,  and  yet  still  able  to  recover 
the  original  signal  x{n)1  It  turns  out  to  be  possible  by  carefully  designing  channel 
filters  [63,  53]. 

For  a  given  frequency  u!k,  STFT  coefficients  of  (2.1)  can  be  rewritten  as 

oc 

X{ejUk,  n)  =  e-iu"n  £  x{m)w(n  -  m)e^("-m)  =  e~iu"n  [x{n)  *  vk(n)} , 

m=— oo 

where  Vk(n)  =  w(n)ejuJkTl.  This  means  that  a  frequency  point  of  STFT  can  be  com- 
puted using  a  filter  with  frequency  response  Vk{uo)  =  W(u  —  Uk)  and  a  complex 
multiplier  e-ju>kU.  Note  that  \X(ejUk,n)\  =  \x(n)  *  vk(n)\. 

For  a  filter  bank  consists  of  M  —  1  such  channels  satisfying 

M-l  M-l 

£  Vk(uj)=  Y  W(u-uk)  =  l,  (2.3) 

k=0  k=0 

the  original  function  x(n)  can  be  recovered  by  replacing  the  integral  of  (2.2)  with  the 
following  summation: 


M-l  M-l 


Y  X(ej"k,n)eju,kn  =  £  [x{n)  *  vk(n)}  =  x(n)  * 


M-l 


Y  vk(n) 


(2.3)     .  . 

=  x{n). 


k=0  k-0  ik=Q 

Figure  2.1  showed  a  filter  bank  structure  for  the  analysis-synthesis  process  of 

STFT. 
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-  vi(n) 


x(n) 


r  "o(n) 


fJW-l(n) 


X"(eJ-o.") 
x)  •  (  x 


A'(e^i-n) 

X)  ■  IX 


x(n) 


Figure  2.1:  A  Filter  Bank  Implementation  of  STFT. 

2.2.6    Generalized  Short-time  Fourier  Transform 

We  already  pointed  out  that  STFT  magnitude  can  be  computed  by  a  filter 
bank  without  modulators  and  demodulators  (see  Figure  2.1).  We  can  further  relax 
the  constraint  that  channel  filters  are  frequency-shifted  of  a  single  lowpass  filter  v(n), 
and  define  a  generalized  short-time  Fourier  transform  pair  as 

00 

wk{n)  —   ^2  x{m)vk{n  —  m)i  0  <    <  M  —  1, 

m=  — oo 
M-l 

x(n)  =  Yl  u'fc(n)  • 

2.3    Wavelet  Transforms 

Wavelet  transforms  have  become  a  powerful  tool  for  analyzing  non-stationary 
signals  in  recent  years.  It  has  several  major  advantages  over  the  short-time  Fourier 
transform. 

First,  wavelet  transforms  choose  to  fix  the  ratio  of  bandwidth/center-frequency 
(called  constant-  Q)  of  channel  frequency  response  instead  of  the  frequency  resolution 
(determined  by  the  bandwidth  of  the  window)  as  in  STFT.  In  other  words,  wavelet 
transforms  use  a  narrow  bandwidth  window  for  low  frequency  signals  and  a  wide 
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bandwidth  window  for  high  frequency  signals.  The  idea  behind  this  is  that  low  fre- 
quency components  represent  slow  changes  in  the  time  domain  for  which  time  resolu- 
tion is  not  critical  while  higher  resolution  in  the  frequency  domain  is  more  desirable. 
In  the  opposite  case,  a  wide  bandwidth  window  for  high  frequency  components  means 
better  resolution  in  the  time  domain  to  capture  abrupt  changes. 

Second,  recent  development  of  wavelet  packet  transforms  and  other  signal- 
adaptive  wavelet-type  transforms  are  much  more  flexible  and  powerful  in  exploiting 
time-frequency  concept  for  wide-ranging  applications. 

2.3.1    Continuous  Wavelet  Transform  (CWT) 

For  a  continuous  time  signal  x(t),  a  continuous  wavelet  transform  pair  may 
be  written  as  [41,  12]: 


wavelet.  Parameter  r  is  the  translation  factor  and  a  is  the  dilation  factor.  Informally 
speaking,  when  a  increases,  function  ipa(t  —  r)  expands  and  takes  long-time  behavior 
into  account.  In  the  opposite,  when  a  decreases,  function  ipa(t  —  r)  contracts  and 
only  focuses  on  the  short-time  behavior.  In  order  to  recover  the  original  signal  from 
its  CWT,  the  wavelet  ip{t)  must  satisfy  the  admissibility  condition: 


However,  a  wavelet  function  tl>a{t)  is  generally  not  an  eigenfunction  of  a  con- 
volution operator  T,  as: 


where  x*  denotes  complex  conjugate  of  x,  ipa(t) 


V>(a)-  ^)  *s  caued  the  basic 


does  not  equal  to  Xtpa(t). 
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Since  inverse  CWT  of  (2.5)  can  be  rewritten  as: 

/•O=0O  (Iq 

x(t)  =  ?r  /       [X(a,t)*i>a(t)}-  , 
Ja=0  cr 

where  *  denotes  convolution  on  variable  t,  we  have: 

1    r°°  da 

Tf[x(t)}  =  —  [(f(t)*X(a,t))*rpa(t)]-f 
O^,  Jo  tr 

If  we  denote  W  as  the  operator  of  continuous  wavelet  transform,  the  above 
result  showed: 

/oo 
f(S)X(a,T -£)<%  =  TfW[x(t)]. 
-oo 

Therefore,  CWT  is  a  translation  invariant  representation. 

2.3.2    Discrete  Wavelet  Transform  (DWT) 

Discrete  wavelet  transform  pair  on  real  orthogonal  wavelet  bases  may  be  writ- 
ten as  [63]: 

/oo 
x{i)2-jl2xl>(2-n  -  n)dt  {DWT), 
-oo 
+00  +oo 

x{t)  =  Y.  E  Vj{n)2->/2ip(2-Jt-n)    {Inverse  DWT). 

j=0  n=-oo 

where  functions  ^2~^2ip{2~H  —  n)|  consist  an  orthonomal  basis. 

Notice  that  DWT  maps  a  continuous-time  function  x{t)  into  a  discrete-time 
sequence  /ij(n).  Clearly,  it  is  not  a  translation  invariant  representation. 

As  in  the  case  of  CWT,  basis  functions  2~^2^{2~H  —  n)  generally  is  not 


a  eigenfunction  of  convolution  operators.  However,  Ty  2  ^2ip{2  H  —  n) 
decomposed  in  DWT  space: 


mav  be 


T, 


2-^(2~H  -  n)\  =  Y,  £  m)2-,/V(2"'t  -  m), 

(=0  m=— oo 


where  the  transform  coefficients  may  be  calculated  by: 

/oo  r  . 

Tf  2-]l24){2~H  -  n)\  2-ll2i>{2~h  -  m)dt 
-oo         L  * 

=  [°°  IT  /(02"J/2^  (2~j{t  -  0  -  n)  del  2-//V(2-'i  -  m)dt 

J-oo  U— oo  '  J 

j -co       yj —oo    v  y 
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x(n) 


9( 


rio(n) 


-®-f$(n7 


-|/i(n)  |— (T2)-  c0(n) 


L|fc(n)[-(|2)-        i  -@-[^) 


cow 


i(n) 


Figure  2.2:  Fast  Discrete-time  Wavelet  Transform-First  Two  Levels. 
Note:  Down  arrow  =  decimation  and  up  arrow  =  expansion. 

We  can  derive  DWT  representation  of  Tj  [x(t)\  to  be 
+00  +00 

Tf[x(t)]  =  E  E  H(n)Tf  [2-"V(2-'«  -  n) 

j=0  n=-oo 


+00  +00 


=  I  E 

/=0  m=— 00 


+c»  +00 


E  E  N(n)viAn,™) 

j=0  n=— oo 


Therefore,  {^(n,  m)}  can  be  seen  as  representation  of  the  operator  Tf  under 
the  basis  {2-^2^{2~H,  -  n)}. 

2.3.3    Discrete-time  Wavelet  Transform  (DTWT) 

The  discrete-time  wavelet  transform  (DTWT)  cannot  be  obtained  by  sam- 
pling the  continuous-time  wavelet  transforms  nor  by  simply  mimicking  the  continuous 
ones.  Assuming  a  discrete  "mother  wavelet"  u(n),  a  discrete-time  version  of  dilation 
v(2kn  —  m)  does  not  produce  a  frequency  response  ~  V(e2  u).  Moreover,  a  dilation 
V{e2ku))  in  frequency  domain  is  generally  not  a  band-pass  function  but  a  multi-band 
function.  The  fundamental  difference  lies  on  both  the  27r-periodic  characteristic  of 
frequency  response  and  the  indexing  constraint  of  discrete-time  series. 

A  fast  algorithm  based  on  the  concept  of  multiresolution  was  found  by  S. 
Mallat  [42]  and  is  illustrated  in  Figure  2.2. 

If  we  consider  only  one  level,  we  have 


Mn)  =  Em=_oo  x(m)g(2n  -  m), 


(2.6) 
(2.7) 
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x(n)  =  ££_oo  co(l)h(n  -  2/)  +  ££_«,  MW"  ~  2/). 


(2.8) 


Therefore,  the  decomposition  filters  h(n)  and  g(n)  and  the  reconstruction 
filters  h(n)  and  g(n)  must  satisfy: 

Em  #(2«  -  m)g(m  -  21)  =6(n-l), 
Em  h{2n  -  m)h(m  -  21)  =  5(n  -  I), 

Zm9(2n-m)h{m-2l)=0,  (2.9) 
£m  h(2n  -  m)g(m  -  21)  =  0, 

£,  /i(2/  -  m)ft(n  -  2/)  +  £,  (7(2/  -  m)g{n  -  21)  =  5(m  -  n) 
It  can  be  proved  (see  Appendix  A)  that  in  frequency  domain  those  conditions 

are  equivalent  to: 

Gie^Gi^)  +  G(^u+1())G{Su+ir))  =  2 
H(e?w)H(e?u)  +  H(ej^+n))H(ej{u}+7[))  =  2 
G{e?u)H{j")  +  G{ej{uJ+ir))H(e^+w))  =  0 
#(e^)G(eJu;)  +  H{e?("+*))G(e>(u+*))  =  0 
H{^)H[eju)  +  G{ejuJ)G(eJul)  =  2 
H(ejul)H(ej{uJ+w))  +  G(eJUJ)G(ej{uJ+ir))  =  0. 

One  possible  choice  is  to  choose  h(n),  g(n),  h(n)  and  g(n)  by: 

g(n)  =  {-l)nh{l  -  n)         G(e>'w)  =  -e~^H*{e^+^), 
h(n)  =  h*(-n)  H(e?»)  =  i/*(eJW), 

g(n)=g*(-n)  #  G(^)  =  G*(e^l 

In  this  case,  the  resulting  bases  is  called  orthogonal  bases.  Otherwise,  the  resulting 


(2.10) 


(2.11) 


bases  is  called  biorthogonal  bases. 

It  can  be  proved  (see  Appendix  A)  that  the  structure  of  Figure  2.2  is  equivalent 
to  a  filter  bank  with  the  transform  written  as  [63]: 

dk(n)  =   ]T  x(m)vk{2k+1n  -  m),   0  <  k  <  M  - 

m=— 00 
00 

CM-i(n)  =   Y.  x{m)vM-\{2M~ln  -  m), 


m=— 00 
M-2  00 


►  (DTWT) 


x(n)  =         X!  dk(m)uk(n  -  2k+lm)  +   ^  cA/-i(m)uA/-i(n  -  2M  1m) 

fc=0  m=— 00  m=— 00 

(/river se  D7WT). 
Particularly,  for  orthonomal  basis,  it  can  be  proved  (see  Appendix  A)  that 

uk(n)  =  v*k(-n), 
uk{m  —  nkl)u*c(m  —  ncn)  —  S(k  —  c)S(l  —  n).  (2-12) 


l.J 


x(n) 


 1  /^~\  Do(n) 

-|   vo(n)    |~(t"oJ  '  (Wj- 


«o(i)  — | 


ui(n)  - 


!>*( 


^)_-@ — — (K)-_M 


(n) 


x(n) 


Figure  2.3:  Generalized  Discrete-time  Filter  Bank  Transform. 

Since  DTWT  can  be  included  into  a  general  framework  to  be  presented  next, 
we  leave  more  in  depth  discussions  to  the  coming  section. 

2.4    Generalized  Discrete-time  Filter  Bank  Transforms 

Previous  sections  showed  that  both  STFT  and  DTWT  can  be  implemented 
by  a  filter  bank  structure.  If  we  ignore  the  modulator  (e~julkri)  and  the  demodulator 
(e^h11)  of  STFT,  the  discrete-time  version  of  the  two  can  be  included  in  a  generalized 
form: 

Mn)  =  E~=-oo  x{m)vk{nkn  -  m),  0  <  k  <  M  -  1,  (GDTFBT)  (2.13) 
x(n)  =  Efclo1  E^=-oc  Mm)uk(n  -  nkm)  {Inverse  GDTFBT).  (2.14) 

It  was  called  the  generalized  filter  bank  transform  [63].  Figure  2.3  illustrated 
the  structure  of  GDTFBT. 

Depending  on  nk,  vk(n)  and  uk(n),  GDTFBT  includes  the  following  special 

cases. 

1)  Orthonomal  wavelet  representations.  nk  =  2fc  for  0  <  k  <  M  —  2  and  Um_i  = 
um-2-  In  this  case,  basis  functions  possess  orthonomality  of  (2.12). 
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-@-fo>(n) 


i(n) 


Figure  2.4:  Example  of  General  Binary-tree-structured  Filter  Bank  of  Five  Channels. 

2)  Orthonomal  wavelet  packet  representations.  This  is  a  generalization  on  the  or- 
thonomal  wavelet  by  applying  the  recursive  decomposition-reconstruction  struc- 
ture to  all  branches.  Therefore,      =  2L  for  0  <  k  <  2L  —  1. 

3)  Biorthogonal  wavelet  (packet)  representations.  The  same  as  the  orthonomal 
counterpart  except  that  the  constraints  (2.11)  are  lifted. 

4)  Generalized  binary  tree  structured  filter  bank  representations.  This  is  a  further 
generalization  upon  the  (bi)orthogonal  wavelet  (packets).  Every  channel  may 
be  split  into  two  with  filters  satisfy  (2.9).  Therefore,  every  is  power  of  2. 
Moreover,  filters  used  in  one  channel  (node)  may  be  different  with  another. 
Figure  2.4  showed  an  example.  It  was  proved  in  Appendix  A  that  the  basis 
functions  constructed  this  way  satisfy  the  biorthonomality  expressed  by: 


uc(m  —  ncl)vki;nf.n  —  m)  =  5(c  —  k)S(l  —  n). 


(2.15) 


5)  Overcomplete  wavelet  representations.  In  this  case,  no  down-sampling  and  up- 
sampling  are  used,  which  means  =  1  for  all  channels  in  the  Figure  2.3.  As 
an  example,  the  overcomplete  wavelet  transform  corresponding  to  Figure  2.4  is 
shown  in  Figure  2.5. 
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Figure  2.5:  The  Overcomplete  Wavelet  Transform  for  the  Example  of  Figure  2.4. 

Representations  1-4  above  are  called  "complete.  The  reason  of  calling  the 
representations  of  the  category  5  "overcomplete"  or  "redundant"  is  that  they  con- 
tain more  analyzing  functions  ({ujt(n)})  than  a  basis.  Since  these  analyzing  and 
synthesizing  functions  must  satisfy 

A/-1      00  A/-1 

£    £  uk(m-l)vk(n-m)=8(l-n),   or    £  Uk{e?«)Vk{e?")  =  1, 

A;=0  m=— 00  k=0 

they  must  be  linear  dependent.  In  the  other  word,  some  functions  must  be  expressible 
by  others. 

Complete  and  overcomplete  representations  behave  rather  differently  under  a 
convolution  operator. 

2.4.1    Convolution  Operator  on  Complete  Representations 

In  general,  a  basis  function  Uk{n— n^m)  is  not  an  eigenfunction  of  a  convolution 
operator.  Using  the  same  idea  as  in  the  case  of  DWT,  a  convolution  operator  Tj  can 
be  represented  by  basis  functions  {u^n  —  n^rn)}: 

M-l  00 

Tf  [uk{n  -  nkm))  =  J2   £  %,c(m>  b)uc(n  -  ncb)  , 

c=0  b=  -00 

00  00 

VkAm' b)  =  £    £  /(')***(*  ~  *  ~  nkTn)vc(ncb  -  I).  (2.16) 

/=— 00  i=— cx> 

It  can  be  proved  that  ujt(n  —  is  an  eigenfunction  of  the  operator  Tj  if 
and  only  if   /7/t,c(m,  6)  =  XS(k  -  c)5{m  —  b). 
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Proof:  If  T]k,c(m^  b)  =  \8(k  —  c)5(m  —  b),  then 

Tj  [Uk(n  -  nkm)}   =   EcL~0l  EZ-oo  ^(k  -  c)S(m  -  b)uc(n  -  ncb) 
=   \uk(n-  nkm). 

In  the  other  direction,  if  Tj  [uk(n  —  nkm)\  =  Xuk(n  —  nkm),  then 

VkAmM   =  En=-ooTf[uk(n-nkm)]vc(ncb-n) 

=   En=-oc  ^uk(n  -  nkm)vc(ncb  -  n)  =  X5(k  -  c)8(m  -  b)  , 

where  we  utilized  the  property  (2.15).  ■ 

The  set  of  coefficients  {rjktC(m,b)}  is  the  representation  of  the  operator  Tj  in 
the  basis  and  is  independent  of  input  signal.  Since  a  convolution  operation  on  a 
sequence  x(n)  can  be  derived  as: 

A/-1  oo 

Tf[x{n)}   =    Y    Y  wk(m)Tf[uk(n-nkm)] 

k=0  m=-oo 

A/-1     oo  M-l  oo 

=    Y    Y  wk(m)  Y   Y  %,c(TO>  b)uc(n  -  ncb) 

k—0  m=  — oo  c—0  b=—oo 


M-l  oo 


=  E  E 

c=0  b=— oo 


M-l  oo 


53    Y  wk(m)riktC(m,b) 


uc{n  -  ncb)  , 


.  k=0  m=-oo 

the  representation  of  T  [x(n)]  is  thus 

A/-1  oo 

acib)=Y    Y  Wk(m)r}k<c(m,b).  (2.17) 

jfc=0  m=-oo 

For  iV-periodic  sequences  x(n),  uk(n)  and  vk(n),  r}k,c(m,b)  is  ^—periodic  in 
m  index  and  ^-periodic  in  b  index,  and  ac(b)  is  ^-periodic  in  b  index. 
A  representation  in  matrix  form  was  derived  by  Banham  et  al  [3]. 
For  a  translation  operator  of  f(i)  =  5(i  —  rf), 

%jC(m, b)  =  E~-oo uk{l  ~d-  nkm)vc(ncb  -  I), 

and  generally  is  not  equal  to  X5(k  —  c)6(m  —  b).  Therefore,  complete  representations 
are  not  translation  invariant. 


1!) 

2.4.2  Convolution  Operator  on  Overcomplete  Representations 

The  overcomplete  representation  of  a  convolution  operator  Tj  can  be  obtained 
from  (2.16)  with  nc  =  1  and  nk  =  1: 

oo  oo 

ilk,c{m,b)=         Y,  f{i)Ml-i-m)vc(b-l)  =  {f*Uk*vc)(b-m).  (2.18) 

/=— oo  :=  —  oo 

Accordingly,  the  overcomplete  representation  of  7/  [x(n)]  can  be  obtained  from 

(2.17): 

M-l  oc 

«c(n)=S    Y  wk{m)r]k,c(n-7n)  =  f(n)*wc{n),  (2.19) 

k—O  m=-oo 

which  showed  that  discrete  overcomplete  wavelet  transform  operator  and  convolu- 
tion operator  are  commutable,  and  thus  overcomplete  wavelet  representations  are 
translation-invariant.  To  verify  this,  assume  an  overcomplete  wavelet  transform  op- 
erator of  W  and  a  translation  operator  of  Tj  with  f(n)  =  S(n  -  d),  overcomplete 
wavelet  representation  of  T)  [x(n)]  is: 

WTj  [x{n)\  =  TfW[x(n)}  =  Tf  [{tufc(n)|o<*<Af_i}]  =  {wk(n  -  rf)|o<*<M-i}, 

which  means  that  the  representation  of  a  time-shifted  signal  is  simply  the  translated 
version  of  the  original  representation  with  the  same  amount  of  time-shift. 

2.4.3  Approximation  of  Convolution  Operators  on  Overcomplete  Repre- 
sentations 

In  frequency  domain,  (2.18)  can  be  written  as: 

Tjk>c(e>»)  =  F{e?")Uk{e>«)Vc{e?")  , 

where  %,c(eJu;)  is  the  discrete-time  Fourier  transform  of  i]k,c(.n)- 

If  the  passbands  of  the  filters  C4(e-7'w)  and  Vc(e^w)  have  little  overlap,  the 
following  approximation  is  justified: 

Uk(enVc(en  «  5(k  -  c)Uc(e^w)Vc(ePu)  . 
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Moreover,  if  the  function  F(eJW)  is  real  and  the  passbands  of  channel  c's  are 
narrow  enough,  we  may  further  approximate  fjk^ie^)  as: 

VkAen  «  AC5(A:  -  c)^e(c>tf)Vc(e>w), 

or, 

%,c(«)  ~  Ac<$(fc  ~  c)uc{n)  *  vc(n)  , 
where  Ac  is  a  real  number.  In  this  case,  (2.19)  may  be  well  approximated  by: 

M-l 

<*c(n)  ~       KHk  -  c)wk{n)  *  uc(n)  *  vc(n)  =  Actuc(n)  .  (2.20) 


k=0 


This  means  that  a  convolution  operator  with  real  frequency  response  F(e*w) 
may  be  approximated  by  a  set  of  real  multipliers  {A^}.  This  is  another  advantage  of 
overcomplete  representation  since  we  found  that  the  approximation  is  generally  not 
extendible  to  complete  representations. 

2.4.4    Aliasing  Enhancement  on  Complete  Representations 

For  complete  representations,  subsampling  usually  does  not  meet  the  Shannon 
sampling  theorem,  and  thus  produces  aliasing  distortions.  Those  aliasing  components 
get  canceled  only  by  specially  designed  filters  at  the  reconstruction  stage.  If  we 
manipulate  the  representation  in  a  way  like  (2.20),  aliasing  distortions  may  not  be 
canceled.  For  a  concrete  example,  we  consider  a  wavelet  transform  of  Figure  2.2 
with  only  one  level.  If  we  multiply  representations  d0(n)  and  c0(n)  by  kd  and  kc 
respectively,  we  can  rewrite  Equations  (2.6-2.8)  in  frequency  domain: 


D'0(e?u)  =  hd  [X(e^/2)G(ej^2)  +  X{e^+2*)/2)G(ej^+2*)/2)}  . 
Cj(e*")  =  ^kc[X{e^'2)H{e^l2)  +  X{e^+2^)l2)H{ej{M)l2)\  , 
X'{e>u)   =   l-  [kdG{ejuJ)G(ejuJ)  +  keH{e^)H(eiu)]  X{e>w)  + 

i  [kdQ{eP{u+v))G(eiw)  +  kcH{ej{uJ+7r))H{ejiJ)]  X(ej(w+,r)). 


The  second  term  is  the  aliasing  component.  Unless  kd  =  kc,  the  second  term 
is  generally  not  equal  to  zero  even  though  filters  meet  the  perfect  reconstruction 
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conditions  (2.10).  In  this  case,  the  system  is  not  a  linear  time  invariant  (LTI) 
system,  but  rather  a  linear  periodically  time  varying  (LPTV)  system  [63].  It  is 
easy  to  recognize  that  the  system  can  not  be  characterized  by  the  familiar  form  of 
Y{e^)  =  M{e^UJ)X{e^ul).  Figure  2.6  shows  an  example  of  such  aliasing-enhancement 
effect. 


50  100  150  200  250  0  0  5  1  1  5  2  2.5  3 

(a)  (b) 


50  100  1  50  200  250  0  0  5  1  1-5  2  2.5  3 


(c)  (d) 

Figure  2.6:  Example  of  the  Aliasing-enhancement  Anomaly. 
Where: 

(a)  The  original  sinusoid  signal; 

(b)  Magnitude  of  DFT  of  the  signal  (a); 

(c)  Reconstructed   signal   with   orthogonal   Haar  filters 

{H{e?u)  =  2-'/2(l  +  e-'u),  G(e>'w)  =  2"1/2(1  -  e"^)), 
one  level  decomposition  and  kj  =  6,  kc  =  1; 

(d)  Magnitude  of  DFT  of  the  signal  (c). 
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2.5    Summary  and  Discussion 

We  have  reviewed  signal  representations  based  on  Fourier  and  wavelet  trans- 
formations and  included  them  in  a  general  framework  called  generalized  filter  bank 
transform.  We  studied  their  behavior  under  a  convolution  operator.  It  revealed  that 
complex  sinusoid  functions  of  Fourier  transform  are  the  only  basis  functions  in  this 
class  which  are  eigenfunctions  of  the  operator.  In  this  case,  a  convolution  operator  is 
equivalent  to  multiplications  in  the  representation  (frequency)  domain  and  the  inverse 
operator  can  be  simply  expressed  as  divisions.  For  overcomplete  representations,  a 
convolution  operator  may  be  approximated  by  multiplications  of  coefficients  by  real 
factors.  However,  such  an  approximation  can  not  be  extended  to  complete  repre- 
sentations due  to  aliasing  distortions.  Moreover,  we  realized  that  discrete  complete 
representations  also  lack  translation-invariant.  Based  on  the  characteristics  of  these 
representations,  we  may  find  the  best  applications  for  them. 

1)  Complete  representations  may  be  good  for  compression,  progressive  transmis- 
sion. No  redundancy  alone  would  be  a  major  advantage  for  those  applications. 
Throwing  away  some  coefficients  may  be  seen  as  a  zeroing  operation,  which 
does  not  enhance  aliasing  components.  After  all,  we  are  prepared  to  accept 
some  distortions  for  those  applications. 

2)  Overcomplete  representations  possess  many  desirable  properties  for  restoration, 
enhancement,  feature  extraction  and  time  delay  estimation  [52].  First,  they  are 
translation  invariant  which  is  critical  for  pattern  recognition.  Second,  they 
do  not  have  the  aliasing  distortions.  Third,  they  allow  much  higher  degree  of 
flexibility  in  filter  design. 


CHAPTER  3 

OVERCOMPLETE  WAVELET  REPRESENTATIONS 
3.1  Introduction 

The  paradigm  of  overcomplete  wavelet  representations  is  a  versatile  and  pow- 
erful tool  shown  to  be  advantageous  for  certain  applications.  In  this  chapter,  we  shall 
study  overcomplete  wavelet  representations  in  detail,  including  issues  of  decomposi- 
tion and  reconstruction  algorithm,  filter  selection  and  time-frequency  interpretation. 

3.2    Overcomplete  Wavelet  Transforms  and  Filters 

In  general,  a  discrete-time  overcomplete  wavelet  transform  (DTOWT)  pair 
can  be  obtained  by  setting  nk  to  1  in  the  Equations  (2.13-2.14),  as  rewritten  in  the 
following: 

wk(n)  =  i>*(ro)a:(n  -  m),  0  <  k  <  M-  1,  (DTOWT)  (3.1) 

x(n)  =  Efcio1  E~=-oo  Mm)uk(n  -  m)  {Inverse  DTOWT).  (3.2) 

More  importantly,  we  need  an  algorithm  to  generate  analyzing  and  synthe- 
sizing filters  Vk(n)  and  uk(n).  In  this  thesis,  we  restrict  ourselves  to  the  binary  tree 
algorithm  as  illustrated  in  the  previous  chapter,  and  consider  two  particular  classes: 
overcomplete  wavelet  packets  and  dyadic  wavelets. 

3.2.1    Overcomplete  Wavelet  Packet  Representations 

The  construction  algorithm  of  overcomplete  wavelet  packet  representations 
corresponds  to  a  complete  binary  tree,  as  shown  in  Figure  3.1.  Though  each  tree 
node  represents  an  analyzing  filter,  only  a  subset  of  all  tree  nodes  is  needed  to  achieve 
a  perfect  reconstruction.  Notice  that  the  position  and  the  width  of  each  rectangular 
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H{2uj)i//'\^G{2uj)  ///\h(2uj)  D 


H{2ui) 

JJ 


#2,0 

1*2,1 

N-2,2 

#2,3 

G{2u)) 


2tt 


(a)  (6) 

Figure  3.1:  Binary  Tree  for  Overcomplete  Wavelet  Packet  Decomposition: 
(a)  the  tree;  (b)  illustration  of  frequency  scaling  of  the  prototype  filters. 

node  also  represents  the  channel's  position  and  the  band  width  in  the  frequency 
domain. 

In  general,  prototype  filters  H(to)  and  G(u)  in  each  level  can  be  distinct,  as 
proved  in  Appendix  A,  as  long  as  they  satisfy  the  condition  H(lu)H(cu)+G(uj)G(lo)  =  1. 
The  frequency  response  of  each  node  can  be  calculated  recursively: 

P0(u)  =  H{u>),  P1(u)  =  G{u)] 

jV0loM  =  1;  Nl+l,m(u)  =  Pc@u)NlAm/21(u),  I  >  0.  (3.3) 

where  the  subscript  c  of  the  filter  Pc(u)  follows  the  periodic  pattern  of  0, 1, 1, 0,  and 
can  be  calculated  by: 

c  =  (m  +  [m/2\)  mod  2. 

The  purpose  of  such  ordering  is  to  make  the  index  m  of  the  node  N^m  proportional 
to  its  central  frequency. 

The  reconstruction  from  two  child  nodes  to  their  parent  node  can  be  expressed 
by  the  following  formula: 

where  we  used  the  property  of: 

Pkmod2{^>)  Pkmod2(^)  +  P(k+l)mod2(u)  P(k+l)mod2(u) 

=    H{lo)H{u)  +  G{lo)G{oj)  =  1. 
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There  are  many  possible  choices  of  filters  Vk{n)  from  the  tree  to  achieve  the  per- 
fect reconstruction.  For  example,  for  the  tree  shown  in  the  Figure  3.1  with  the  height 
of  two  [55,  page  449],  five  sets  are  available:  {A^o},  {Ni,o,  iV"i,i},  {^1,0,^2,2,^2,3}, 
{^2,0,^2,1,^1,1}  and  {N2,o,  #2,1'  N2t2,N2,3}.  In  general,  let  T(h)  be  the  number  of 
selections  from  a  complete  binary  decomposition  tree  of  height  h.  Since  the  tree  may 
be  split  into  two  subtrees  of  height  h-l  each  in  turn  have  T(h-l)  selections,  we 
have  the  following  recursive  formula: 

T(h)  =  [T{h  -  l)]2  +  1,  h>l,   and  T(0)  =  1. 

The  first  five  numbers  were  calculated  as  T(l)=2,  T(2)=5,  T(3)=26,  T(4)=677  and 
T(5)=458330. 

All  quadrature  mirror  filters  (QMF),  good  for  orthogonal  wavelet  transform, 
are  also  good  for  overcomplete  wavelet  packet  transforms.  However,  for  real  QMF 
except  Haar  wavelet,  I.  Daubechies  proved  that  linear  phase  and  finite  impulse  re- 
sponse (FIR)  properties  are  mutual  exclusive  [13].  For  example,  Daubechies  wavelet 
filters  [11]  are  FIR  but  not  linear  phase.  In  the  other  hand.  Battle-Lemarie  wavelet 
filters  are  symmetric  but  not  FIR,  which  may  be  written  as  [42]: 


Hp(u) 

where  p  is  a  positive  integer  and 


(3.4) 


+00  j 

Note:  An  algorithm  for  numerical  computation  of  Battle-Lemarie  wavelets  is  included 
in  Appendix  B. 

3.2.2    Dyadic  Wavelet  Representations 

The  building  structure  of  a  dyadic  wavelet  transform  is  a  particular  subset 
of  the  overcomplete  wavelet  packet  transform.  It  decomposes  low-frequency  branch 
only  (assume  H{uj)  is  lowpass),  as  shown  in  Figure  3.2. 
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Figure  3.2:  Binary  Tree  for  Dyadic  Wavelet  Decomposition. 

Obviously,  all  filters  good  for  overcomplete  wavelet  packet  transforms  are  good 
for  dyadic  wavelet  transforms.  However,  one  may  want  to  exploit  the  extra  oppor- 
tunity created  by  the  more  relaxed  constraint  on  filters.  The  following  filter  classes 
[44,  20]  possess  distinct  property  in  the  time  domain: 


H(u)  =  ejpu/2cos2n+p{uj/2), 
H{J)  =  H*(u), 

G(u)  =  -(-j)Pe-^/2sin2-p(u;/2), 

G{u)  =  -UY^smnul2)Y^Tlca^{ul2) 


(3.5) 


where  n  >  0  is  an  integer  and  p  £  {0, 1}. 

In  fact,  for:  p  =  0  filter  g(-l)  =  0.25,  g(0)  =  -0.5  and  g(l)  =  0.25  is  propor- 
tional to  a  discrete  Laplacian  operator  and  for:  p  =  1,  g(0)  =  0.5  and  #(1)  =  -0.5  is 
proportional  to  a  discrete  gradient  operator. 

Notice  that  filters  within  this  class  are  linear  phase  and  FIR,  but  not  orthog- 
onal. For  p  =  0,  filters  H,  H,  G  and  G  and  their  frequency  scalings  {e.g.  H(2lu)) 
are  all  symmetric.  For  p  =  1,  different  prototype  filters  may  be  used  for  frequency 
scaling  to  minimize  spatial  shifting  caused  by  the  phase  factor.  If  we  denote  Hi(u), 
Hi(u>),  G[(lu)  and  Gi(to)  for  filters  used  in  level  /  >  1,  we  find  the  following  filters  are 
either  symmetric  or  antisymmetric  as  well  as  FIR: 

Ht(u)  =  cos2"+1(2'-2o;),   Hfa)  =  Ht{u), 

2n 

Gt(u)  =  j  sin(2'-2u;),       Gi(u)  =  -Gt{u)  £  cos2m(2'-2o;). 


m=0 
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Figure  3.3:  Examples  of  the  Function  Qa,b(u)- 

In  both  cases,  frequency  response  of  tree  node  N^k  (I  >  0  and  k  G  {0, 1})  can 
be  derived  as: 


l-i 


Nlfi{u)  =  e^'2  J]  cos2n+p(2m-1u;)  =  e 


jpu/2 


Nu{u)  =  -(-j)pe 


771  =  0 


2'-1sin(u;/2) 


2-P 


sin(2'-1u;) 
2'  sin(w/2) 
sin(2'-2o;) 


2n+p 


2n+2 


2'"1  sin(u;/2) 

For  convenience,  we  defined  a  function  Qa<b(u))  as: 


©o,ft(w) 


sin(2Q-1tj) 


2a  sin(w/2)_ 

For  6  >  4,  the  function  QaA00)  ls  approximately  Gaussian  as  plotted  in  Figure  3.3. 
Similarly,  we  can  define  a  reconstruction  tree  with  nodes: 

Nifi(u)  =  iV/_i,o(w)G,(w). 
A  dyadic  wavelet  representation  consists  of  leaf  nodes.  For  example,  a  three- 
level  representation  consists  of  nodes  {N3fi,  jV3)1,  AT2)i,  Ni,i}.  In  terms  of  filter  bank, 
they  can  be  viewed  as  multiple  channels.  Channel  frequency  responses  can  be  derived 

as: 

Clfi(u))     =     Nlfi{iO)Nlfi{uj)  =   Yl  Hm(u)Hm(u)  =  9l,4n+2p(u), 

771=1 

CM(w)   =   JVu(w)iV,il(a;)  =  iV/_1,oMNi_1,o(u;)G,(a;)G/(a;) 
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=   iV^oMA^oM  [l  -  H^H^uj) 

=     @i-l,47i+2p(^')  —  ©/,4n+2p(^)- 

The  channel  frequency  responses  for  a  three-level  representation  are 
{C3fl,C3,i,C2,i,C\%i}.  Notice  that  the  channel  frequency  responses  Ctti(u)  are  an 
approximate  difference  of  two  Gaussians  (DOG)  [47,  page  63].  Two  examples  of 
channel  filters  are  shown  in  Figure  3.6. 

3.2.3    Filters  of  Autocorrelation  Functions  of  Wavelets 

From  a  filtering  point  of  view,  the  idea  is  very  simple.  Since  there  is  no 
downsampling-upsampling  taking  place  in  overcomplete  wavelet  transforms,  it  seems 
that  we  can  simply  combine  decomposition  and  reconstruction  into  a  single  filtering 
step.  To  achieve  this,  one  only  need  to  change  prototype  filters  into: 

Ha{u)  =  H{u)H(u),  (36) 
Ga(co)  =  G(u)G{u). 

The  perfect  reconstruction  condition  now  takes  the  form  of: 

Ha(u)  +  Ga(u)  =  l. 

As  H(u)  is  connected  to  a  wavelet,  Ha{uj)  is  connected  to  an  autocorrelation 
function  of  wavelet.  For  rigorous  mathematical  treatment,  refer  to  Saito's  Ph.D  thesis 
[57]. 

In  this  case,  a  reconstruction  step  is  simply  a  summation  of  all  coefficients  as 
(3.2)  degenerates  to: 

x(n)  =  ZkJo1  ufc(n), 

which  takes  the  same  form  as  the  generalized  STFT. 

The  most  important  property  is  that  filters  Ha(u)  and  Ga(uj)  are  always  sym- 
metric or  antisymmetric,  and  thus  providing  a  simple  way  to  construct  symmetry 
and  FIR  filters  for  certain  applications. 


29 


3.3    Connection  to  Complete  Wavelet  Representations 

In  Appendix  A,  Proposition  A. 0.3  built  the  bridge  between  an  overcomplete 
wavelet  representation  and  a  complete  wavelet  representation.  If  both  representa- 
tions used  the  exactly  same  prototype  filters  H(u),  G(u),  H(uj)  and  G(u),  the  only 
difference  lies  on  the  decimation.  In  this  case,  the  complete  wavelet  representation 
is  a  subset  of  the  corresponding  overcomplete  wavelet  representation  due  to  decima- 
tion. If  the  prototype  filters  used  by  an  overcomplete  wavelet  representation  does  not 
satisfy  (2.10),  there  is  no  corresponding  complete  wavelet  representation. 

3.4    Time-frequency  Interpretation  and  the  Uncertainty  Principle 

The  set  of  transform  coefficients  {wk(n)\0<k<M_l}  is  a  representation  of  the 
signal  x(n)  in  the  transform  domain.  It  is  important  to  point  out  that  an  analyzing 
filter  Vk{co)  with  nonzero  phase  response  will  cause  a  time  shift  [49,  page  205]  dk  in 
the  corresponding  component  Wk{n).  For  some  applications,  e.g.,  segmentation,  this 
shift  has  to  be  compensated.  In  this  case,  {wk(n  +  <4)|0<fc<M-i}  should  be  used  as 
the  representation. 

Figure  3.4  shows  two  overcomplete  wavelet  packet  representations  using  sym- 
metric Lemarie  filter  of  p  =  1.  It  is  apparent  that  overcomplete  wavelet  representa- 
tions include  both  frequency  and  temporal  information  of  the  signals.  The  frequency 
resolution  or  capability  to  separate  different  frequency  component  of  an  overcom- 
plete wavelet  representation  is  determined  by  Vfe(u;)'s.  The  narrower  the  passband 
of  the  Vk(u)  is,  the  higher  it  is  the  frequency  resolution.  On  the  other  hand,  the 
time  resolution  is  determined  by  the  distribution  (or  window  size)  of  the  impulse  re- 
sponses i>fc(n)'s.  The  shorter  the  distribution  is,  the  higher  the  time  resolution  it  can 
achieve.  Unfortunately,  high  resolution  in  time  and  frequency  domains  are  conflict 
goals,  as  demonstrated  in  the  Figure  3.5.  In  this  example,  we  saw  that  as  pass- 
band  narrowed,  impulse  responses  spread  and  consequently  the  boundary  between 
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Figure  3.4:  Two  Examples  of  the  Overcomplete  Wavelet  Representation. 
Legend: 

(a)  An  ideal  square  waveform; 

(b)  A  linear  chirp  signal  s(n)  =  5cos(0.027rn  +  0.00l7rn2). 

In  each  image,  top  half  is  the  signal  waveform,  and  the  bottom 
half  is  (0-255)  scaled  overcomplete  wavelet  representation  (coef- 
ficients) of  64  leaf  nodes  of  the  level  six.  For  the  bottom  half, 
horizontal  axis  is  time,  vertical  axis  is  frequency. 
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two  impulses  blurred.  This  phenomenon  reflected  a  fundamental  principle  called  the 
uncertainty  principle  [21,  51]. 

Time  and  frequency  localization  of  a  discrete  lowpass  function  f(n)  can  be 
mathematically  described  by  two  quantities, 

^  =  4£(n-n)2!/(n)|2, 

n 

where  E  =  £„  |/(n)|2  =  £  £,  |F(e*)|a<k;  ,  and  n  =  J  £n  "|/(«)|2- 
The  product  of  <r^  and  er2  is  defined  as  the  uncertainty  factor. 

U  =  a2nol  (3.7) 

Liu  and  Akansu  [40]  proved  that  U  >  [E  -  \F(ejn)f]  /(4£2),  and  thus  for  filters  with 
F(ejn)  =  0,  U  >  0.25.  That  the  uncertainty  factor  has  a  lower  bound  is  significant 
and  is  called  the  uncertainty  principle.  It  basically  states  that  we  cannot  have  filters 
with  arbitrarily  narrow  bandwidth  in  the  frequency  domain  and  arbitrarily  short 
duration  in  the  time  domain.  As  a  result,  we  cannot  achieve  arbitrarily  precise  time 
and  frequency  localization  simultaneously.  (Numerical  computation  of  uncertainty 
factors  for  filter  banks  is  discussed  in  Appendix  C.) 

Figure  3.6  shows  examples  of  analyzing  filters  of  overcomplete  wavelet  trans- 
form and  their  uncertainty  factors.  Uncertainty  factors  of  the  analyzing  filters  of  the 
dyadic  wavelet  showed  that  they  are  indeed  very  close  to  Gaussian  functions.  In  gen- 
eral, the  results  showed  that  uncertainty  factor  varied  from  channel  to  channel,  and 
level  to  level.  As  a  statistical  measure,  we  used  maximum,  minimum  and  mean  to 
describe  higher  levels.  Such  statistics  for  higher  level  overcomplete  wavelet  packets 
are  presented  in  the  Table  3.1,  where  Regular  stands  for  regular  Lemarie  filters  of 
(3.4)  and  Autocorrelation  stands  for  autocorrelation  functions  of  Lemarie  filters. 

We  may  catch  some  trends  from  the  data  on  Lemarie-filter-spanned  analyzing 
filters  of  overcomplete  wavelet  packet  representations: 
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Col  1  Col  2  Col  3 


Figure  3.5:   Examples  of  Time-frequency  Localization  by  Overcomplete  Wavelet 

Packet  Representations. 

Where: 

Row  1:  A  signal  consists  of  two  pure-tone  impulses  with  frequency  of  u  =  0.3n 
and  u!  =  0.327r:  the  spectral  (Col  1)  and  the  waveform  (Col  3). 

Row  2-5:  level  2  representation. 

Row  6-13:  level  3  representation. 

For  row  2-13:  columnl,  frequency  responses  {A^ifc|0  <  k  <  2'};  column2,  impulse 
responses  of  columnl  (dot  lines  approximate  envelopes).  column3,  overcomplete 
wavelet  representation.  The  prototype  filter  used  was  Lemare  filter  of  p  =  1. 
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Row  1 


Col  1 


Col  2 


Figure  3.6:  Uncertainty  factors  (drawn  on  top  of  each  band)  of  Channel  Filters. 
Where: 

Row  1:  Overcomplete  wavelet  packet  of  level  2  using  Lemarie  filter  of  p  =  1  (Col 
1)  andp  =  9  (Col  2); 

Row  2:  Overcomplete  wavelet  packet  of  level  3  using  Lemarie  filter  of  p  =  1  (Col 
1)  and  p  =  9  (Col  2); 

Row  3:  Channel  filters  {C4j0,  C4)i,  C3)i,  €2,1,  C^i}  of  dyadic  wavelet  using  the  fil- 
ter class  of  (3.5)  with  n  =  1  (Col  1)  and  p  =  1,  n  =  1  (Col  2). 
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Table  3.1:  Uncertainty  Factors  of  Analyzing  Filters 


Regular 

Autocorrelation 

p 

level 

T  J 

"max 

J  j 

T  T 

U 

^max 

C-'min 

TT 
L 

4 

1 .618 

0.350 

0.701 

0.8o3 

c\  or™ 

0.36  ( 

0.49< 

1 

1 

0 

O.Ool 

U.04Z 

i  7At; 
1. 1  uo 

z.uoo 

U.OIO 

6 

16.721 

0.318 

4.535 

5.693 

0.262 

0.980 

4 

1.011 

0.437 

0.628 

1.178 

0.470 

0.584 

2 

5 

3.464 

0.397 

0.954 

0.871 

0.333 

0.432 

6 

7.984 

0.320 

1.809 

2.736 

0.260 

0.475 

4 

1.283 

0.529 

0.674 

1.349 

0.496 

0.646 

3 

5 

2.386 

0.418 

0.704 

0.721 

0.329 

0.398 

6 

5.244 

0.318 

1.005 

1.605 

0.261 

0.346 

1)  Joint  time- frequency  localization  of  those  using  low  order  filters  tends  to  dete- 
riorate quickly  as  level  increase.  This  is  due  to  the  existence  of  small  sidelobes. 

2)  Using  filters  of  autocorrelation  function  is  more  effective  than  using  higher  order 
regular  filter,  if  the  no-reconstruction  structure  is  acceptable  by  applications. 

However,  one  should  be  very  careful  in  using  uncertainty  factors  as  an  optimality 
measure  due  to  two  main  reasons.  First,  how  can  we  come  up  a  single  sensible 
number  for  a  filter  bank?  Minimum  Umax  or  minimum  U  or  something  else?  Second, 
more  fundamentally,  the  uncertainty  factor  does  not  take  the  signal  into  account.  For 
configurations  with  similar  uncertainty  factors,  as  in  the  case  of  Gaussian-like  filters 
such  as  (3.5),  the  measure  of  uncertainty  factor  alone  would  not  tell  us  which  one  is 
optimal. 

3.5    Two-dimensional  Extensions 

So  far  we  have  only  discussed  one-dimensional  representations.  In  order  to 
apply  overcomplete  wavelet  transforms  to  two-dimensional  images,  we  need  to  extend 
the  transform  to  two-dimension.  The  following  two  separable  extensions  can  be 
implemented  as  separated  ID  filtering  of  rows  and  columns. 


35 


1)  Tensor  product  extension  [42].  This  extension  is  originally  applied  to  orthogonal 
wavelet  transform.  Only  ID  filters  H(u)  and  G(u>)  are  used.  It  directly  exploit 
the  perfect  reconstruction  property  of  quadrature  mirror  filters  by  extending: 

H{u))H{u)  +  G{u)G{u>)  =  1 

into: 

H(ux)H{ux)  +  G(ux)G{ux)]  [H{uy)H(Luy)  +  G(ujy)G(ujy) 

The  above  2D  extension  can  be  rewritten  as 

HH(ujx,  ujy)HH(ujx,LUy)  +  HG{ux,  uy)HG(u)x,u}y) 

+ 

GH(ux,ujy)GH(ujx,ujy)  +  GG{ujx,u)y)GG(u)x,u}y)  =  1. 

where 

HH(ujXlUy)  =  H{ux)H{ujy),  HH{ujx,ujy)  =  H(u)x)H{ujy), 
HG{lux,  ujy)  =  H(ujx)G(ujy),  HG(u;x,  ujy)  =  H(ujx)G{uy), 
GH{ux,uy)  =  G(ujx)H(ujy),  GH{jjx,ujy)  =  G(ux)H(ujy): 
GG(ujx,uy)  =  G(ux)G{ujy),  GG(ujx,ujy)  =  G(ujx)G{uiy). 

Filters    {HH,  HG,  GH,  GG}    are    used    for    decomposition    while  filters 
{HH,  HG,  GH,  GG}  are  used  for  reconstruction.  The  same  construction  method 
is  applied  to  frequency  scaling  of  higher  levels  with  each  node  decomposed  into 
four  child  nodes. 

2)  Two-orientational  extension.  For  the  particular  dyadic  wavelet  filters  (3.5)  with 
emphasis  on  spatial  operation  (gradient  and  Laplacian),  the  component  {GG} 
of  the  tensor  product  method  may  not  be  useful.  Therefore,  the  following  2D 
extension  was  proposed  by  Mallat  [44]: 

HH(ujx,ujy)  =  H{ujx)H{ujy),    HH(cux,u!y)  =  H(UX)H(iUy), 
GX(ux,ujy)  =  G(lux),  GX(u)x,uy)  =  G(uJx)L{uy), 

GY(ux,ujy)  =  G(uy),  GY(vx,vy)  =  L(ux)G(uy). 

This  2D  extension  has  only  three  components: 


36 


•  a  horizontal  component  GX 

•  a  vertical  component  GY 

•  a  lower  resolution  component  HH 

However,  it  introduced  a  new  filter  L(u),  which  is  to  be  determined  by  the 
condition  of: 


HH  HH  +  GX  ■  GX  +  GY  ■  GY  =  1. 


It  turned  out  that  the  formula  for  L(uj)  is 


L{u) 


l  +  H(u;)H{u) 
2 


For  filters  of  (3.5), 


L(u;)  =  0.5  *  1  +  cos 


,4n+2p 


(w/2)]  , 


which  is  a  symmetric  and  FIR  filter. 


CHAPTER  4 

OPERATIONS  ON  OVERCOMPLETE  WAVELET  REPRESENTATIONS 

4.1  Introduction 

In  most  cases,  a  suitable  signal  representation  is  not  the  end  of  a  problem- 
solving  process.  On  the  contrary,  how  to  manipulate  it  to  achieve  a  desired  effect  is 
the  more  challenging  part. 

For  overcomplete  wavelet  representations,  we  have  identified  three  operations: 
gain,  shrinkage  and  envelope  operators  for  our  applications  in  image  deblurring  and 
segmentation.  We  demonstrated  that  a  linear  convolution  may  be  approximated  by 
a  gain  operator,  and  devised  a  tree  pruning  algorithm  based  on  a  tolerant  of  approxi- 
mation error.  We  extended  the  wavelet  shrinkage  proposed  by  Donoho  and  Johnstone 
[17]  to  overcomplete  wavelet  representations  for  denoising.  We  built  the  connection 
between  the  power  density  distribution  of  a  overcomplete  wavelet  representation  and 
its  envelopes,  and  presented  two  envelope  detection  algorithms. 

4.2    Gain  Operators 

In  the  Section  2.4.3,  we  determined  that  a  convolution  operator  can  be  approx- 
imated by  multiplication  of  real  factors  {Afc|0<jt<M-i},  followed  by  a  reconstruction. 
Such  a  operation  can  be  denoted  by  a  gain  operator  Q.  When  applied  to  a  over- 
complete  wavelet  representation,  a  gain  operator  Q  multiplies  Wk(n)  by  a  real  gain 
of  Afc.  Notice  that  such  gain  operators  are  linear  and  memoryless.  We  called  the  set 
of  gain  factors  { Afc |o<fc<A/-i}  £  $M  a  gain  vector.  In  fact,  such  a  operation  with  an 
arbitrary  gain  vector  is  equivalent  to  a  linear  filtering.  For  a  given  channel  configura- 
tion and  a  gain  vector  {A^  |o<fc<A/— i  }i  the  relationship  between  a  signal  x(n)  and  the 
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reconstructed  signal  x'(n)  in  the  frequency  domain  is  simply 

M-1  M-1 

X'(u)  =  Y,  hWk{u>)Uk(u)  =  X(u)  Y  \kVk(u)Uk(u), 

k=0  k=0 

and  thus  the  equivalent  filter  is 

Q(w)  =  [X0Xi  ■  ■  •  AM-i]  [C0{u)d(u)  ■  ■  •  CM-,(w)f  , 

where  [A0Ai  •  •  •  Am-i]  is  a  gain  vector,  Ck(u;)  =  Vk(u)Uk(u)  is  the  channel  filters,  and 
Ml  denotes  the  matrix  transpose  of  M. 

Two  fundamental  issues  immediately  arise: 

1)  What  is  the  characteristics  of  the  set  of  filters  expressible  by  a  gain  operator?  In 
the  other  words,  what  kind  of  filters  may  be  approximated  by  a  gain  operator? 
A  broad  answer  is  that  any  gain  operator  can  only  represent  a  symmetric  filter. 
This  is  because  channel  filters  {Ck(uj)\0<k<M -1}  of  overcomplete  wavelet 
representations  are  usually  real  functions,  any  linear  combination  of  them  must 
also  be  a  real  function. 

2)  For  a  given  frequency  response  F(u),  how  can  we  determine  a  tree  configuration 
and  a  gain  vector  as  the  best  approximation?  Note  that  this  question  really 
asks  for  a  criteria  of  "best"  and  an  algorithm  to  find  it. 

4.2.1    Minimum  Mean  Square  Error  (MMSE)  Approximation 

The  closeness  of  the  approximation  in  terms  of  frequency  response  by  a  wavelet 
representation  with  channel  frequency  responses  {Cfc(u;)|o<jt<A/-i}  and  a  gain  vector 
{Afc|o<fc<M-i}  to  a  given  real  function  F(u)  may  be  measured  by  the  following  mean 
square  error. 


1  r 

"  2?r  J-* 


M-1 


k=0 


duo.  (4.1) 
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It  is  well  known  that  the  optimal  gain  vector  corresponding  to  the  minimum 
mean  square  error  may  be  found  by  setting  all  partial  derivatives       to  zero: 


35         1  r 


1  r* 

7T  J-TT 


M-l 


k=0 


Cq{uo)duj  =  0,  0  <  q  <  M-l, 


which  is  a  set  of  linear  equations: 

M-l 


Y  IT  Ck{u))Cq(u)duj\  \k  =  r  F{Lu)Cq(u>)duu,  0  <  q  <M-1.  (4.2) 


1       rir  M~!         1  rn 

Srmn  =  ^  /     F*(u>)du  -  £  A£—  /  F '{u)Ck{u>)<b> , 
Z7T  J-tt  ,_n        zTT  J -w 


k=0 

The  minimum  mean  square  error  can  be  derived  as 

«r  M    1  pit 

k=0 

where  {A£|o<9<a/-i}  is  the  optimal  gain  vector  satisfying  (4.2). 

We  may  define  a  normalized  minimum  mean  square  error  (NMMSE)  e  as 

6min 

where  Ej  =  j-  f*n  F2(u)du  is  the  energy  of  the  given  frequency  response.  Since 
&min  >  0,  it  must  be  true  that 

Ef  >  E^Lo1       £r  F(u)Ck(u)du, 

and  therefore  0  <  e  <  1  . 

The  criteria  of  mean  square  error  provides  a  way  to  determine  the  optimal 
gain  vector  for  a  given  frequency  response  F(u)  and  a  given  configuration.  However, 
the  criteria  alone  suggests  that  the  bottom-most  nodes  are  the  best  representation, 
for  which  e  may  reach  its  minimum  among  all  possible  configurations.  This  can  be 
formally  proved  by  the  following  proposition. 

Proposition  4.2.1  Splitting  a  channel  into  two  will  not  increase  the  error  measure 
e. 
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Proof:  By  contradiction.  Assume  that  for  the  current  configuration  {Cfc(w)|0<fc<A/-i}, 
the  optimal  gain  vector  is  {X°k}  and  the  NMMSE  is  e.  We  further  assume  that  the 
channel  q  is  split  into  two  channels  with  indexes  qx  and  q2,  with  the  optimal  gain 
vector  {X°k}  for  the  new  configuration  and  the  new  e'  >  e. 
Without  loss  of  generality,  we  can  write: 

Cqi(u)=Cq(u)H(2lu,)H(2lu>), 

CQ2(u)  =  Cq{uj)G{2lu;)G{2luj).  K  > 

Since  the  prototype  filters  satisfy  H(lo)H(u))  +  G(w)G(uj)  =  1,  we  have 
Cqi(u))  +  Cq2{u)  =  Cq(uj).  For  the  new  configuration,  we  can  find  a  gain  vector 
{\'k  =  \°k,  k  /  qi,  k^q2\  \qi  =  X'q2  =  A°}  with  error  e'  =  e.  This  contradicts  the  assump- 
tion that  e1  >  e,  and  therefore  e'  <  e  ■ 

Examples  showing  higher  level  corresponding  to  better  approximation  can  be 
visualized  in  Figure  4.1.  Overcomplete  wavelet  packet  representations  of  leaf  nodes  at 
Level  4,  5,  6,  7  were  selected,  corresponding  to  16,  32,  64,  128  channels,  respectively. 
Lemarie  filter  of  p  =  3  were  used. 

4.2.2    Time-frequency  Trade-off  and  a  Greedy  Algorithm 

Proposition  4.2.1  clearly  states  that  approximation  error  e  of  frequency  re- 
sponse decreases  as  the  channel  bandwidth  decreases,  as  demonstrated  by  the  exam- 
ples in  Figure  4.1.  Therefore,  by  the  criteria  alone,  the  Fourier  transform  would  be 
the  best  representation  since  its  approximation  error  is  zero.  This  is  certainly  not 
the  answer  we  are  looking  for. 

What  we  lost  in  sight  is  the  other  part  of  the  picture  in  the  time  domain. 
As  dictated  by  the  uncertainty  principle,  time  resolution  decreased  as  the  frequency 
resolution  increased.  A  overcomplete  wavelet  packet  representation  by  the  deepest 
leaves  would  be  closed  to  a  Fourier  representation  and  thus  lost  its  resolution  in 
the  time  domain.  Moreover,  from  computational  point  of  view,  such  a  down-to-the- 
bottom  decomposition  is  rather  inefficient. 
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Figure  4.1:  Examples  of  Gain  Vector  Approximation  (frequency  range  [0,  it]  shown) 
Where: 

Row  1:  a  given  frequency  response  F(eJ1J); 

Row  2  to  5:  approximated  frequency  responses  by  overcomplete  wavelet  packet 
representations  of  level  4,  5,  6,  7,  overlaid  with  F(e^)  (dotted  curve),  with 
5=39.4%,  18.5%,  4.5%,  0.02%,  respectively. 
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One  possible  trade-off  between  time  and  frequency  can  be  accomplished  by 
setting  up  a  tolerant  on  the  approximation  error  e.  Under  an  error  tolerant,  the 
goal  could  be  to  find  a  tree  configuration  with  minimum  number  of  leaves.  In  the 
following,  we  present  a  greedy  algorithm  in  the  order  of  breadth-first  search  [45]. 
Notice  that  a  tree  node  here  represents  a  channel  filter  C^m(u)  and  the  relationship 
between  children  nodes  and  their  parent  node  is  characterized  by  (4.3). 

Algorithm  4.2.2  (Minimum  Tree)  Given  a  frequency  response  F(u>)  and  an  er- 
ror tolerant  £q,  let  current  leaves  C  =  {Co,o(w)  =  1}  and  compute  the  current  error 
£c  using  the  set  C. 

Step  1.  One  by  one,  compute  error  Ek  by  temporarily  substituting  a  node 
in  the  current  leaves  C  with  its  two  children.  Choose  the  pair  with  minimum  Ek  to 
replace  their  parent  in  the  set  £  and  update  the  current  error  Ec  =  m'm{Ek}- 

Step  2.  If  Ec  <  Eq,  stop  and  use  the  leaves  L  as  representation;  else,  go  to 
step  1. 

Examples  of  overcomplete  wavelet  packet  approximation  using  the  minimum 
tree  algorithm  are  shown  in  Figure  4.2.  The  improvement  is  appreciated  by  compar- 
ing it  with  Figure  4.1.  In  Figure  4.1,  the  representation  using  thirty-two  channels 
of  the  Level  5  has  the  approximation  error  s  —  18.5%.  However,  for  the  same  fre- 
quency response,  the  minimum  tree  algorithm  was  able  to  find  a  representation  of 
only  twenty-six  channels  and  yet  reduced  the  approximation  error  to  e  =  6%.  In  the 
other  case,  the  minimum  tree  algorithm  found  a  representation  with  merely  seven 
channels  and  achieved  approximation  error  of  0.6%. 

4.3    Shrinking  Operators 

Wavelet  shrinkage  was  introduced  by  Donoho  and  Johnstone  [17]  for  signal 
denoising  and  linear  inverse  problems.  For  a  thorough  introduction,  refer  to  [15]. 
The  description  of  the  methodology  was  given  in  [15,  page  173],  we  quote,  "Wavelet 
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Row  1 


Col  1 


Col  2 


Figure  4.2:  Examples  of  the  Minimum  Tree  Approximation  with  the  Maximum  Depth 

of  Seven  Levels. 

Where: 

Col  1:  A  given  frequency  response  Fi(ejuJ)  (rowl,  same  as  in  Figure  4.1),  ap- 
proximation result  overlaid  with  Fx  (row2)  with  e  =  6%  and  selected  26 
channels  (row3); 

Col  2:  A  given  frequency  response  F2(eju')  (rowl),  approximation  result  overlaid 
with  F2  (row2)  with  e  =  0.6%  and  selected  7  channels  (row3). 

Overcomplete  wavelet  packet  with  Lemarie  filter  Hi(ejuJ)  was  used  in  both  cases,  and 
the  frequency  range  [0, 7r]  is  shown. 
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shrinkage  refers  to  reconstructions  obtained  by  wavelet  transformation,  followed  by 
shrinking  the  empirical  wavelet  coefficients  towards  zero,  followed  by  inverse  trans- 
formation." For  a  thorough  introduction,  refer  to  [15]. 

We  extended  the  idea  to  overcomplete  wavelet  representations  on  the  ground 
that  they  differs  with  orthogonal  wavelet  representations  mainly  on  having  extra 
correlated  coefficients. 

A  shrinking  operator  K.t  can  be  written  as 


Kt  [w] 


w  —  t,    w  >  t, 

0,  \w\  <  t,  (4.4) 

w  +t,    w  <  —t, 


where  t  >  0  is  a  threshold.  Notice  that  the  shrinking  operator  is  a  nonlinear  and 
memoryless  operator,  and  is  supposed  to  work  in  spatial  domain. 

The  goal  of  a  shrinking  operation  is  to  cut  the  noise  off  while  keeping  useful 
signals.  Obviously,  this  is  possible  only  if  the  amplitude  of  desirable  components  is 
stronger  than  the  noise  components.  Therefore,  the  threshold  value  t  is  dependent  on 
the  signal-to-noise  ratio.  Examples  of  choosing  appropriate  thresholds  are  discussed 
in  Chapter  6. 

To  demonstrate  the  shrinking  operation  and  the  advantage  of  wavelet  denois- 
ing,  an  example  is  included  in  Figure  4.3.  The  major  advantage  of  wavelet  shrinkage 
denoising,  as  compared  with  linear  smoothing,  is  that  it  diminishes  noise  but  does  not 
smear  the  sharp  edges,  as  seen  in  the  example.  Comparison  with  a  linear  averaging 
filter  with  h  =  [1, 1, 1, 1, 1, 1,  l]/7  is  presented  in  Figure  4.4. 

4.4    Envelope  Detectors 

Energy  distribution  is  a  very  important  concept  in  the  physical  world.  It  is  a 
indispensable  tool  for  signal  analysis,  too. 

For  a  harmonic  current  i(t)  =  As'm((jj0t  +  a)  passing  through  a  one-ohm  pure 
resistant,  the  instantaneous  power  (energy  per  unit  time)  at  time  t  is  i2(t).  The 
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Col  1 


Col  2 


Figure  4.3:  An  Example  of  Overcomplete  Wavelet  Denoising. 
Where: 

Col  1:  An  perfect  impulse  contaminated  by  white  Gaussian 
noise  (rowl)  and  its  level  2  overcomplete  wavelet 
packet  representation  (row2-5); 

Col  2:  Overcomplete  wavelet  packet  representation  after  shrink- 
ing operations  (row2-5)  and  reconstructed  signal 
(rowl). 
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Row  1  •* 
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2  :: 


Figure  4.4:  Comparison  of  Overcomplete  Wavelet  Shrinkage  Denoising  with  Linear 

Smoothing  Using  the  Signal  of  (row  l,col  1)  of  Figure  4.3. 

Where: 

Row  1:  Overcomplete  wavelet  shrinkage  (left,  the  same  as  row  1,  col  2 
of  Figure  4.3)  and  linear  smoothing  (right); 

Row  2:  Local  zoom-in  with  overlay  of  the  original  noisy  signal. 

average  normalized  power  in  a  period  T  is  thus 

W=\[i\t)dt=\A\ 

Therefore,  the  average  normalized  power  of  a  harmonic  signal  is  proportional 
to  its  amplitude  scjuared. 

For  an  arbitrary  waveform  i(t),  one  can  decompose  it  into  harmonic  compo- 
nents using  Fourier  analysis: 

/oo 
I(u)e>2*ftdf. 
-oo 

By  the  same  token,  energy  density  (per  unit  time  and  per  unit  frequency)  at  frequency 
uj  is  proportional  to  the  amplitude  squared  of  the  component  7(a;)eJ'27r^t,  which  is 
|/(u;)|2.  The  function  |/(<x>)|2  describes  power  distribution  in  the  frequency  domain 
and  is  called  power  density  spectrum  [49,  page  66]. 

The  above  concepts  can  be  extended  to  overcomplete  wavelet  representations, 
by  extending  the  concept  of  amplitude  to  "envelope"  [26,  chapter  4],  or  "local  ampli- 
tude." An  overcomplete  wavelet  representation  is  consist  of  bandpass  components. 
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Since  any  bandpass  waveform  may  be  represented  by  A(t)  sin(ucf  +  a(t))  [26,  page 
229],  where  A(t)  is  the  envelope  and  uic  is  the  associated  carrier  frequency,  an  over- 
complete  wavelet  representation  {wjfc(t)|o<fc<A/-i}  maY  be  rewritten  in  the  form  of: 


Thus,  squared  envelopes  {A2(uk,t)}  shall  describe  the  power  density  distribution 
of  the  overcomplete  wavelet  representation  in  the  time-frequency  plane  (also  called 
"phase  space"  in  [12]).  Figure  4.5  shows  examples  of  such  a  distribution. 

For  time-frequency  representations  using  complex  Gabor  filters  [6],  envelopes 
may  be  extracted  by  simply  performing  a  modulus  operation  on  two  quadratic  com- 
ponents. However,  for  representations  with  real  analyzing  filters,  such  as  overcom- 
plete wavelet  representations,  more  sophisticated  envelope  detection  algorithms  are 
needed.  We  present  the  next  two  envelope  detection  algorithms  and  investigate  their 
performance. 

4.4.1    Envelope  Detection  by  Hilbert  Transform 

The  envelope  of  a  narrowband  bandpass  signal  can  be  computed  by  a  cor- 
responding analytical  signal  [51].   For  a  signal  x(t),  the  analytic  signal  is  defined 


The  envelope  of  the  original  signal  x(t)  is  then  simply  the  modulus  of  the 
analytic  signal  x(t): 


The  frequency  characteristics  of  the  Hilbert  transform  can  be  expressed  by: 


{A(uk,t)sm(ukt  +  ak(t))\o<k<M-i}  ■ 


by: 


x{t)  =  x(t)+jx{t), 


where  x(t)  is  the  Hilbert  transform  of  x(t), 


e(t)  =  \x{t)\  =  y/x*(t)+S?(t). 


uj>=0 
otherwise. 
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(b) 

Figure  4.5:  Two  Examples  of  Power  Density  Distribution  of  Overcomplete  Wavelet 

Packet  Representations. 

Legend: 

(a)  A  linear  chirp  signal,  the  same  as  in  Figure  3.4(b); 

(b)  A  multi-segment  pure  tone  signal  with  frequencies  uj0  =  0.087T,  u>i  =  0.257T, 
o>2  =  0.157T  and  u3  =  0.557T,  in  the  order  of  occurrence.  The  overcomplete 
wavelet  representation  and  visual  arrangements  are  the  same  as  in  Figure  3.4. 
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Figure  4.6:  A  FIR  Hilbert  Transformer. 


Legend: 

(a)  the  imaginary  part  of  the  frequency  response  (the  real  part  is  zero). 

(b)  nonzero  coefficients. 

Therefore,  the  Fourier  transform  of  the  analytical  signal  x(t)  is: 


For  implementations  in  discrete-time  domain,  approximate  FIR  Hilbert  trans- 
formers may  be  designed  by  windowing  the  ideal  frequency  response  [49].  Figure  4.6 
shows  an  example.  It  is  a  type-Ill  FIR  Hilbert  transformer  designed  with  parameters 
M  =  18,  and  (3  —  2.629  [49.  page  680].  This  Hilbert  transformer  is  antisymmetric,  of 
length  19  and  has  only  10  nonzero  coefficients. 

4.4.2    Envelope  Detection  by  Zero  Crossings 

In  this  method,  the  maximum  absolute  value  between  two  adjacent  zero- 
crossings  is  first  found,  and  then  assigned  to  each  point  within  the  interval. 

Algorithm  4.4.1  (Envelope  by  Zero  Crossings)  For  a  given  array x(l  :  N),  start 
from  the  index  i  =  1,  find  the  next  index  k  which  is  either  a  zero  crossing  point  or 
a  zero  valued  point  x(k)  =  0  or  the  other  end  point  k  =  N,  whoever  is  first  encoun- 
tered, assign  the  maximum  absolute  value       =  max  {|x(n)|j<„<fc}  to  all  the  elements 


uj>=0 
otherwise. 
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{x(n)\i<n<k} ■  Advance  the  index  i  =  k  and  repeat  the  process  until  reaching  the  other 
end. 

4.4.3  Comparison  Between  the  Two  Detectors 

We  compared  two  envelope  detectors  by  their  working  frequency  range  and 
robustness  under  noise  perturbations.  For  pure-tone  impulse  signals,  both  did  well  in 
the  middle  range  (u>  €  [0.17T,  0.77r]).  While  Hilbert  transform  was  unacceptable  at  the 
low  end  (u  <  O.ln),  the  zero-crossing  approach  did  poor  at  the  high  end  (uj  >  0.77r). 
The  test  cases  were  shown  in  the  rows  1  and  2  of  Figure  4.7.  We  constructed  a  noisy 
impulse  by  introducing  a  random  frequency  perturbation: 

s{n)  =  sin((w0  +  0.05  *  RANDN(n))n), 

where  RANDN  (a  MATLAB  function)  is  a  random  noise  with  normal  distribution, 
and  filtered  it  with  a  Lemarie  filter  Hi(u).  This  case  was  included  in  the  row  3 
of  Figure  4.7.  We  found  that  the  zero-crossing  method  exhibited  edge-preserving 
smoothing  characteristics  and  was  more  robust  to  wide-band  noise. 

4.4.4  Comparison  with  Other  Energy  Operators 

A  useful  and  simple  energy  operator  was  analyzed  on  [46].  The  operator  \& 
for  a  discrete  signal  x(n)  was  given  as 

*  [x(n)\  =  x2(n)  -  x(n  +  l)x(n  -  1). 

However,  this  operator  differs  with  envelope  detectors  in  the  way  that  its  value  on  a 
pure  tone  signal  is  not  only  proportional  to  the  amplitude  squared  but  also  a  function 
of  frequency  [46] : 

*  [A  sin(u0n  +  a)}  =  A2  sin2(u;o). 

4.4.5  Two  Dimensional  Envelope  Detection 

Next  we  extended  the  ID  envelope  detection  algorithms  for  the  analysis  of 
two-dimensional  image  signals.  In  the  frequency  domain,  a  two-dimensional  analytic 
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Col  1  Col  2  Col  3 

Figure  4.7:  Comparison  of  the  Two  Envelope  Detectors. 

Where: 

Col  1:  Input  signals  (top  row,  u0  —  0.1ir,  Wi  =  0.57r;  center  row, 
a;o  =  0.057r,  a;i=0.757r;  bottom  row,  noisy  signal  with 
a;0  =  0.l7r). 

Col  2:  Envelopes  detected  via  the  19-tap  Hilbert  transformer. 
Col  3:  Envelopes  detected  via  the  zero-crossing  method. 


signal  may  be  obtained  by  setting  an  appropriate  half  plane  to  zero,  based  on  its 
orientation.  That  is,  for  a  2D  signal  f(x,y),  the  Fourier  transform  of  an  analytic 
signal  f(x,y)  is  either: 


F(uJX,UJy) 


_  f  2F(u)x,uy)    ,   ujx  >-  0 
I  0  ,  otherwise 


or, 


F(u>x,uy)  =  l  2F^y)  >  u»>=° 


0  ,  otherwise. 

For  the  2D  filters  constructed  by  the  tensor  product  method  with  halfband  filters, 
the  equivalent  complex  quadrature  filters  exhibited  the  frequency  response  shown  in 
Figure  4.8.  Notice  that  lowpass  and  diagonal  components  may  have  an  alternate 
arrangement  by  zeroing  out  the  left  or  bottom  half  plane. 

This  separable  property  allowed  one  to  compute  the  envelope  of  a  2D  signal 
using  the  ID  algorithms  described  earlier  in  a  straightforward  manner,  described  in 
the  following. 
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Figure  4.8:  Frequency  Response  of  Equivalent  Complex  Quadrature  Filters  (level  1). 
Legend: 

(a)  lowpass  channel 

(b)  vertical  channel 

(c)  horizontal  channel 

(d)  diagonal  channel 

The  diagonal  shadowed  areas  identify  zeroed  half  planes. 

Algorithm  4.4.2  (2D  Envelope  Detection)  For  a  given  2D  array  w(m,n)  and 
its  orientation.  If  its  orientation  is  horizontal,  apply  the  ID  envelope  detector  column- 
wise; If  its  orientation  is  vertical,  apply  the  ID  envelope  detector  row-wise;  Otherwise, 
apply  the  ID  envelope  detector  column-wise. 


CHAPTER  5 
APPLICATION  I:  TEXTURE  SEGMENTATION 

5.1  Introduction 

In  this  chapter,  we  shall  apply  overcomplete  wavelet  representations  to  the 
problem  of  segmentation  of  textured  images. 

However,  it  seems  that  we  were  back  to  the  beginning.  It  may  sound  odd,  but 
it  is  true  that  there  is  not  even  a  precise  definition  of  "texture,"  let  alone  "texture 
segmentation"  [61].  In  spite  of  the  difficulty,  there  is  a  general  consensus  that  texture 
is  a  property  of  area.  We  felt  that  the  concept  of  texture  is  parallel  to  the  concept  of 
"local  frequency"  for  an  non-stationary  signal.  The  concept  of  local  frequency  is  easy 
to  understand.  It  describes  the  rate  of  change  occurs  in  a  window  around  a  particular 
time.  Yet  there  is  not  a  unique  definition  of  local  frequency.  Therefore,  our  treatment 
of  texture  segmentation  is  to  embrace  it  into  the  paradigm  of  time-frequency  analysis. 

There  have  been  many  research  works  on  texture  segmentation  using  the  con- 
cept of  spatial-frequency  representation.  The  analyzing  functions  used  by  researchers 
include: 

•  Laws  microtexture  masks  (filters)  [24] 

•  Complex  prolate  spheroidal  sequences  [64] 

•  Complex  Gabor  functions  [6,  19] 

•  Real  Gabor  functions  [30] 

•  The  Wigner  distribution  [54] 

More  recently,  wavelet  and  wavelet  packet  representations  have  been  added  to 
the  list  [34,  35,  36,  62,  9,  37]. 

A  segmentation  process  usually  consists  of  two  distinct  phases:  feature  ex- 
traction and  clustering.  Features  for  texture  representation  are  of  crucial  importance 
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for  accomplishing  segmentation  [27].  In  this  chapter,  we  demonstrated  that  over- 
complete  wavelet  packet  representation  and  envelope  detection  make  a  good  feature 
extraction  scheme  on  a  variety  of  both  natural  and  synthetic  textured  images. 

5.2    Texture  Feature  Extraction 

The  texture  features  we  chose  are  the  envelopes  of  overcomplete  wavelet  rep- 
resentations. The  interpretation  of  the  feature  is  that  their  squares  corresponds  to 
spatial-frequency  power  density  distribution  as  pointed  out  in  the  last  chapter.  The 
feature  extraction  thus  consists  of  two  stages:  1)  overcomplete  wavelet  packet  decom- 
position and  2)  envelope  detection. 

For  2D  envelope  detection  on  overcomplete  wavelet  packet  representations,  we 
classified  each  node  in  the  decomposition  tree  into  four  possible  categories,  taking 
into  account  orientation,  as  follows: 

•  The  root  node  is  omnidirectional. 

•  The  node  last  filtered  by  Gi(u>x)Hi(u>y)  corresponds  to  vertical- orientation. 
(Highpass  filter  G<  is  applied  row-wise  and  lowpass  filter  Hi  column-wise.) 

•  The  node  last  filtered  by  Hi(ux)Gi(u)y)  corresponds  to  horizontal- orientation. 
(Lowpass  filter  Hi  is  applied  row- wise  and  highpass  filter  Gi  column-wise.) 

•  The  node  last  filtered  by  Gi(ux)Gi(uy)  corresponds  to  diagonal- orientation. 
(Highpass  filter  Gt  is  applied  row-wise  and  highpass  filter  Gi  column-wise) 

•  The  node  last  filtered  by  Ht(ujx)Ht(uy)  has  the  same  orientation  as  its  parent. 
(Lowpass  filter  Ht  is  applied  row-wise  and  lowpass  filter  Ht  column-wise) 

At  the  end  of  feature  extraction,  each  sample  of  a  signal  was  attached  a  fea- 
ture vector.  In  our  case,  such  a  feature  vector  is  consist  of  envelopes  of  overcomplete 
wavelet  packet  representation.  It  is  possible  to  apply  a  monotonic  function  to  the 
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feature  space  to  transform  it  into  another  feature  space.  For  example,  a  square  func- 
tion x2  will  transform  the  envelope  space  into  a  space  of  power  density  distribution 
and  log  function  log(x)  will  transform  the  envelope  space  into  a  space  of  log  spectral. 


We  raised  the  issue  of  filter  selection  in  the  Chapter  3  on  overcomplete  wavelet 
representations.  For  the  application  of  segmentation,  we  argued  that  symmetry, 
frequency  response,  and  boundary  accuracy  are  important  factors  in  the  selection  of 
filters  for  feature  extraction.  In  the  following,  we  discuss  these  constraints  in  terms 
of  overall  performance. 

•  Symmetry.  For  accomplishing  texture  segmentation,  accuracy  in  texture 
boundary  detection  is  crucial  for  reliable  performance.  In  this  application, 
filters  with  symmetry  or  antisymmetry  are  clearly  favored.  Such  filters  have 
a  linear  phase  response,  where  the  delay  (shift)  is  predictable.  Alternatively, 
filters  with  nonlinear  phase  may  introduce  complex  distortion.  Moreover,  sym- 
metric or  anti-symmetric  filters  are  also  advantageous  in  alleviating  boundary 
effects  through  simple  methods  of  mirror  extension  (see  Appendix  D). 

•  Optimal  Frequency  Response.  In  order  to  derive  an  ideal  filter  frequency 
response  for  the  chosen  feature,  we  considered  a  two-band  filter  bank  with  input 
signals  of  infinite  length  consisting  of  two  segments  with  distinct  pure  tones. 
The  input  signals  can  be  written  as: 


where  .4i>0  and  A2>0.  Except  for  the  boundary  (n  =  0),  we  derived  the 
feature  vectors  (envelopes  of  channel  outputs)  as, 


5.3    Considerations  on  Filter  Selection 


A\  cos{u)in  +  Oil)  ,n<0, 
A2  cos(tu2n  +  012)   ,  n  >  0, 


T  =  (eH,eG) 


fleft  =(A1\H(oj1)\,A1\G(u1)\)  ,n<0, 
fright  =  (A2  \HM\,A2  \G(u2)\)   ,n  >  0. 
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The  angle  9  between  vectors  fieft  and  Tright  is 

\HMH(u2)\  +  \G{an)GM\ 


cos  1 


Vl^OI2  +  \GM\2y/\H(u2)\2  +  \G(u2)f 


and  is  bound  by  O<0<  f .  The  distance  between  the  two  classes  in  the  feature 
space  is  then 


D 


Ti 


left 


+ 


-  right 


■  it  ft 


r, 


right 


cos{9) 


For  9  =  | ,  vectors  fieft  and  fright  are  orthogonal,  and  the  distance  D  reaches 
its  maximum.  Clearly,  the  maximum  distance  in  feature  space  between  any  two 
classes  is  optimal  for  segmentation  and  classification  applications,  in  the  sense 
that  the  classes  are  better  seperated  and  more  robust  to  noise  perturbations. 
Notice  that  for  H{u>x)  +  0  and  G{uj2)  ±  0,  cos(0)  =  0  if  and  only  if  G{u)i)  =  0 
and  H(u2)  =  0.  This  means  that  optimal  filter  banks  should  have  no  overlap 
in  the  frequency  domain,  and  thus  filter  H(u)  with  optimal  frequency  response 
must  be  a  perfect  half-band  filter.  Indeed,  [50]  derived  a  similar  optimal  filter 
for  image  coding  applications. 

Of  course  such  ideal  filters  cannot  be  realized  in  practice.  Based  on  the  preced- 
ing discussion,  filter  bands  overlapping  in  the  frequency  domain  tend  to  reduce 
class  distance  in  feature  space.  Therefore,  a  filter  H(uj)  with  a  large  stop-band 
attenuation  and  flat  pass-band  frequency  response  is  desirable. 

Spatial  Localization.  There  are  two  types  of  boundaries  for  signals  of  finite 
length:  1)  boundaries  of  regions  exhibiting  distinct  characteristics,  and  2)  the 
physical  boundaries  of  a  data  segment.  Feature  vectors  close  to  the  boundaries 
will  be  affected.  The  size  of  the  affected  region  depends  on  the  length  of  the 
channel  filters,  and  the  distribution  of  filter  coefficients.  Therefore,  filters  of 
short  length  and  fast  decay  shall  be  better  suited  for  boundary  detection. 


57 


Unfortunately,  the  preceding  three  criteria  cannot  be  satisfied  simultaneously. 
As  previously  pointed  out,  quadrature  mirror  filters  (QMF)  with  compact  support 
cannot  be  symmetric  or  anti-symmetric.  This  means  that  a  symmetric  constraint  on 
any  QMF  will  be  in  conflict  with  spatial  localization  goals.  Moreover,  large  atten- 
uation in  the  stop  band  requires  a  longer  filter,  which  in  turn  degrades  the  filter's 
spatial  localization.  Again,  we  faced  the  limitation  of  the  uncertainty  principle,  and 
the  previous  discussions  about  the  uncertainty  factor  of  a  filter  bank  is  applicable 
here. 

In  this  study,  we  chose  the  Lemarie  filters  for  its  symmetry  and  good  frequency 
characteristics. 

5.4    The  Basic  Isodata  Clustering  Algorithm 

Segmentation  algorithms  accept  a  set  of  features  as  input  and  output  a  con- 
sistent labeling  for  each  pixel.  Basically,  this  is  a  multi-dimensional  data  clustering 
problem  with  no  general  algorithm  available  [23].  Clustering  algorithms  that  have 
been  previously  used  for  texture  segmentation  can  be  divided  into  two  categories: 
1)  supervised  segmentation  and  2)  unsupervised  segmentation. 

In  practice,  unsupervised  segmentation  is  often  desirable  and  easy  to  vali- 
date. It  is  particularly  useful  in  those  cases  where  testing  samples  are  difficult  to 
prepare  (making  supervised  segmentation  infeasible).  Thus,  we  used  a  Baisc  Isodata 
clustering  algorithm  [18,  page  201]. 

Algorithm  5.4.1  (Basic  Isodata)  Given  a  2D  image  array  x  of  structure  contain- 
ing feature  vectors  and  label  fields,  and  the  number  of  classes  Nc, 

Step  1.  Scan  the  representation  matrix  x  in  raster  order.  For  every  pixel 
encountered,  randomly  pick  a  label  from  set  {0,  ...,NC  —  1}  and  assign  it  to  the  pixel. 

Step  2.  Compute  the  class  center  {C^  |o<A:<ivc-i }  by  calculating  the  mean 
vector  for  each  class  k. 
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Step  3.  Rescan  the  whole  representation  matrix  x,  and  assign  pixel  to 
the  class  k  if  the  Euclidean  distance  between  the  feature  vector  of  the  pixel  and  the 
class  center  Ck  is  the  closest. 

Step  4.  //  no  pixel  changes  its  class  in  the  Step  3,  stop;  else,  go  to  the 
Step  2. 

This  algorithm  differs  from  a  K-means  algorithm  [29]  in  only  one  aspect:  Ba- 
sic Isodata  updates  class  centers  after  a  complete  scan  of  an  input  feature  set  while 
K-means  updates  for  every  reassignment.  In  our  experiments,  Basic  Isodata  outper- 
formed K-means  for  almost  all  cases. 

5.5  Postprocessing 

The  Basic  Isodata  algorithm  labels  each  pixel  independently  and  does  not  take 
into  account  the  high  correlation  between  neighboring  pixels.  A  postprocessing  stage, 
such  as  the  relaxation  labeling  [25],  can  be  used  to  incorporate  some  neighborhood 
constraint  into  the  segmentation  process. 

For  simplicity,  we  used  median  filtering  as  our  postprocessing  to  simulate  the 
benefit  of  a  local  constraint.  In  particular,  a  5x5  median  filtering  was  repeatedly 
applied  to  an  initial  segmentation  until  no  change  of  labeling  occurred. 

5.6    Experimental  Results 

Our  segmentation  algorithm  was  tested  on  both  one-dimensional  signals  and 
two-dimensional  textured  images.  Our  test  images  included  samples  of  two  distinct 
families  of  textures.  In  all  examples  shown,  straightforward  envelopes  of  the  over- 
complete  wavelet  packet  representation  were  used  as  texture  feature. 

5.6.1    One-dimensional  Signals 

Figure  5.1  shows  the  segmentation  of  a  signal  consist  of  a  sinusoid  segment 
and  a  triangular  segment  with  the  same  period  (sixteen  samples)  and  amplitude. 
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Figure  5.1:  Segmentation  of  a  ID  Signal  Consists  of  Triangular  Waveform  and 
Sinusoid. 
Where: 

Row  1:  the  original  signal  and  the  segmentation  result; 

Row  2-5:  the  overcomplete  wavelet  packet  representation  and  envelopes. 


With  envelopes  of  four  channels  of  the  Level  2  representation,  the  perfect  result 
was  achieved.  Figure  5.2  is  a  segmentation  example  of  a  noisy  signal  consist  of  two 
sinusoid  segments  (u  =  0.3n  and  u;  =  0.57r)  contaminated  by  white  Gaussian  noise. 
Envelopes  of  eight  channels  of  the  Level  3  representation  was  able  to  achieve  the 
perfect  result. 

5.6.2    Natural  Textures 

Here  we  used  textures  obtained  from  the  Brodatz  album  [7]  and  public  archive. 
Each  testing  sample  was  first  histogram-equalized  so  that  a  segmentation  result  based 
only  on  first  order  statistics  was  not  possible.  Experimental  results  are  displayed  in 
Figure  5.3.  Experimentally,  we  observed  that  a  lower  order  Lemarie-Battle  filter 
(p=l)  performed  well  in  boundary  detection  (Figure  5.3  (b)),  while  the  higher  order 
Lemarie-Battle  filter  (p=2)  did  a  better  job  within  non-boundary  (internal)  regions. 
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Figure  5.2:  Segmentation  of  a  Noisy  ID  Signal  Consists  of  Two  Pure  Tone  Segments. 
Where: 

Row  1:  the  original  signal  and  the  segmentation  result, 

Row  2-9:  the  overcomplete  wavelet  packet  representation  and  envelopes. 


Col  1 


Col  2 


Col  3 


Figure  5.3:  Segmentation  Results  of  a  Image  Consist  of  Natural  Textures. 
Where: 

Row  1:  Test  image  Tl  (256x256,  8-bit)  consists  of  D17,  herringborn 
weave  and  bark  (true  boundary  overlaid  for  display  only). 

Row  2:  Clustering  results  using  features  extracted  from  a  Level  4 
filter  bank  generated  by  (Col  1)  Lemarie-Battle  filter  of  p  =  1, 
(Col  2)  autocorrelation  function  of  Lemarie-Battle  filter  of 
p=l,  (Col  3)  Lemarie-Battle  filter  of  p  =  2.  The  zero-crossing 
algorithm  were  used  for  envelope  detection. 


Row  3:  Final  segmentations  after  postprocessing. 
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5.6.3    Synthetic  Textures 

We  tested  the  performance  of  our  algorithm  on  several  texture  images  synthe- 
sized from  random  noise  based  on  [27]. 


•  Gaussian  Lowpass  and  Bandpass  Textures.  The  images  were  generated 
by  filtering  a  white  Gaussian  noise  with  mean  of  zero  and  standard  deviation 
of  30  with  a  Gaussian  filter  of  frequency  response  (in  polar  coordinate): 


•  Filtered  Impulse  Noise  (FIN).  The  images  were  generated  by  filtering  a 
random  impulse  image  I(m,  n)  generated  by: 


where  RAN  G  (0.0,1.0)  is  a  random  number  with  uniform  distribution,  with 
Gaussian  impulse  response: 


In  both  cases,  images  were  linearly-scaled  into  the  range  of  [0, 255]. 

Figure  5.4  shows  a  segmentation  result  on  a  Gaussian  lowpass  texture  im- 
age, and  Figure  5.5  shows  a  segmentation  result  on  a  filtered  impulse  noise  (FIN) 
texture  image.  For  this  difficult  test  image  T3,  the  algorithm  achieved  outstanding 
performance. 

We  also  tested  our  algorithm  on  a  texture  image  containing  regular  and  sparse 
elements.  Figure  5.6  demonstrates  the  accurate  segmentation  result. 


A  quantitative  comparison,  presenting  the  accuracy  of  our  segmentation  re- 
sults of  textured  images  is  summarized  in  Table  5.1.  This  performance  is  consistent 


G{r,e)=e 


-(6-e0)2/(2^BZ)e-(r-Fc)2/(2S2r) 


g(m,n)  =  e 


m2/(2S2)e-n2/(2S2)_ 


5.7    Summary  and  Discussion 


(a)  (b)  (c) 

Figure  5.4:  Segmentation  of  a  Synthetic  Gaussian  Lowpassed  Texture  Image. 
Legend: 

(a)  Test  image  T2  (256  x  256,  8-bit):  Gaussian  LP,  left:  isotropic 
Fc  =  0,  Sr  =  60;  right:  non-isotropic  Fc  =  0.  5r  =  60,  90  =  0, 
B0  =  0.175. 

(b)  Initial  segmentation  (Lemarie-Battle  filter  of  p  =  1,  Level  4  and 
zero-crossing  envelope  detection) 

(c)  Final  segmentation  after  postprocessing. 


(a)  (b)  (c) 

Figure  5.5:  Segmentation  of  a  Synthetic  FIN  Texture  Image. 
Legend: 

(a)  Test  image  T3  (256  x  256,  8-bit):  Filtered  impulse  noise,  left: 
non-isotropic  T  =  0.15,  Sx  —  1.0.  Sy  =  1.5;  right:  non-isotropic 
T  =  0.15,  Sx  =  2.0,  Sy  =  1.0. 

(b)  Initial  segmentation  (Lemarie-Battle  filter  of  p  =  1,  Level  4  and 
zero-crossing  envelope  detection) 

(c)  Final  segmentation  after  postprocessing. 


64 


(a)  (b)  (c) 

Figure  5.6:  Segmentation  of  a  Synthetic  Texture  Image  with  Line  Patterns. 
Legend: 

(a)  Test  image  T4  (256  x  256,  8-bit) 

(b)  Final  segmentation  (Lemarie-Battle  filter  of  p  =  1,  Level  3  and 
zero-crossing  based  envelope  detection). 

(c)  Detected  boundary  overlaid  with  the  original  image. 
Table  5.1:  Boundary  Accuracy  of  the  Segmentation  Results. 


Test  image 

Maximum  ABE 

Average  ABE 

±a 

Tl 

11.0 

3.2 

2.6 

T2 

15.0 

2.9 

2.4 

T3 

1  1.0 

2.9 

2.6 

ABE:  Absolute  Boundary  Error  (in  pixels). 


with  the  difficulty  of  segmentation  perceived  by  human  observers.  We  observed  that 
boundary  errors  are  dependent  on  shape  i.e.,  complex  boundaries  yielded  higher 
variance. 

Overcomplete  wavelet  packet  representations  and  envelope  features  provide  a 
versatile  and  flexible  framework.  The  examples  showed  that  satisfactory  results  can 
be  achieved.  Based  on  our  experiments,  low-order  Lemarie  filter  (p  =  1)  achieved  the 
best  results.  However,  we  are  aware  that  there  were  several  important  parameters 
to  be  determined  manually  for  the  segmentation  algorithm  to  work.  These  include 
the  configuration  of  the  representation  and  the  number  of  classes.  On  the  face  of 
lacking  a  precise  definition  of  the  problem,  we  believed  that  no  mathematical  solution 
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is  possible,  and  these  parameters  can  only  be  designed  case-by-case  for  particular 
applications. 

Comparing  with  Gabor  filters,  overcomplete  wavelet  packet  representations 
possess  several  potential  advantages  for  further  exploration: 

•  channel  filters  cover  exactly  the  frequency  domain  (provide  a  mathematically 
complete  representation)  without  significant  overlap,  thus  greatly  reducing  cor- 
relations between  features  extracted  from  distinct  channels 

•  adaptive  pruning  of  a  decomposition  tree  makes  possible  the  reduction  of  com- 
putational complexity  and  the  length  of  feature  vectors 

•  by  moving  up  the  decomposition  tree,  feature  vectors  at  distinct  resolution 
levels  can  be  obtained  to  increase  the  accuracy  of  boundary  detection 


CHAPTER  6 
APPLICATION  II:  IMAGE  DEBLURRING 

6.1  Introduction 

Image  deblurring  is  a  special  case  of  signal  restoration.  The  goal  is  to  recover 
an  image  which  is  blurred  and  usually  degraded  by  noise.  During  the  past  two 
decades,  large  amount  of  research  works  have  been  done  in  this  area  [33,  2,  28]. 

Unlike  the  texture  segmentation,  image  deblurring  is  a  well  defined  problem. 
The  mathematical  model  commonly  used  to  describe  the  degradation  process  is  a 
convolution  and  an  additive  noise,  written  as: 

y(ku  k2)  =  d{ku  k2)  *  x(kuk2)  +  n(kuk2), 

where  x(ki,k2)  is  the  original  image,  y(ki,k2)  is  the  observed  image,  d(ki,k2)  is 
the  impulse  response  of  the  blurring  operator,  n(kx,k2)  is  the  noise  and  the  symbol 
*  stands  for  convolution.  The  deblurring  problem  is  to  restore  the  original  image 
x(ki,k2)  from  the  degraded  image  y(h,  k2).  This  is  thus  a  linear  inverse  problem. 

At  the  first  glance,  the  problem  seems  solvable.  Assuming  known  impulse 
response  of  the  blurring  process  and  the  absent  of  noise,  one  may  formally  derive  the 
solution  in  the  frequency  domain  as: 

X(ui,oj2)  =  Y(ujulo2)/D(uji,uj2). 

This  is  the  simplest,  most  intuitive  approach  and  is  called  inverse  filtering  [2,  28,  22]. 
However,  the  frequency  response  D(u\,u)2)  associated  with  a  blurring  process  usually 
has  zeros  in  the  frequency  plane  and  its  reciprocal  is  thus  unbounded.  In  the  other 
point  of  view,  blurring  generally  causes  irreversible  information  loss  even  without  the 
existence  of  noise.  By  taking  into  account  of  noise,  the  inverse  filtering  approach 
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would  produce  a  restoration: 

X'(ui,u}2)  =  Y(ui,u}2)/D{ui,u2)  =  X(ui,u2)  +  N{u)i,u)2)/D(u)i,u)2)- 

In  this  case,  the  closeness  of  the  restored  image  X'  to  the  original  image  X  is  hinged 
on  the  strength  of  the  noise  N.  Notice  that  the  precise  waveform  of  the  noise  is 
unknown  and  thus  cannot  be  completely  subtracted  from  the  observed  image  Y. 

Inverse  problems  associated  with  a  unbounded  inverse  operator  is  said  to  be 
ill-posed  [33,  14].  The  traditional  mathematical  tool  for  dealing  with  such  problems 
is  regularization  methods  [60]. 

In  this  chapter,  we  shall  apply  overcomplete  wavelet  packet  representations  to 
the  deblurring  problem.  Our  approach  explicitly  suppresses  noise  using  the  nonlinear 
shrink  operator  and  approximates  an  inverse  filter  using  the  gain  operator. 


In  this  section,  we  shall  review  modified  inverse  filters,  Wiener  filters  and  the 
more  recent  theory  of  wavelet-vaguelette  inversion. 

6.2.1    Modified  Inverse  Filters 

For  blurs  having  zeros  in  the  frequency  plane,  there  are  two  modifications 
available  to  make  a  bounded  inverse  filter. 

1.  Gain-limited  Inverse  Filter.  The  idea  is  simply  setting  a  upper-bound  for 
the  magnitude  of  frequency  response  to  grow,  as  described  below: 


where  R  is  a  positive  real  value  and  D(ux,uiy)  is  a  real  filter. 

2.  Pseudo-inverse  Filter[28,  page  276].  This  is  the  same  as  the  gain-limited 
inverse  filter  except  the  constant  gain  region  set  to  zero: 


6.2    Review  of  Some  Deblurring  Techniques 


'  1/R 

F{UX,Uy)    =    <  l/D{LUX,UJy) 


0  <  D(ux,uy)  <  R, 

\D{WX,Uy)\  >  R. 

-R  <  D(ux,uy)  <  0, 


(6.1) 


0 

\/D(UJX,U)y) 


\D(u)x,u)y)\  <  R, 
\D{u)x,uy)\  >  R, 


(6.2) 
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where  R  is  a  positive  real  value. 


Both  approaches  have  their  pros  and  cons.  The  advantage  of  the  gain-limited 
inverse  filter  is  that  it  is  continuous  in  the  frequency  domain.  The  disadvantage 
is  that  the  passband  of  \D(u}x,uy)  \  <  R  may  significantly  enhance  wideband  noise. 
In  opposite,  the  pseudo-inverse  filter  will  completely  wipe  out  anything  outside  its 
region  of  inversion  by  sacrificing  the  continuity. 

6.2.2    Wiener  Filters 

The  idea  of  the  Wiener  filter  is  to  choose  a  deblurring  filter  /  to  minimize  the 
mean  square  restoration  error  defined  as  [18,  2,  22]: 

£  =  E{[x'(kuk2)-x(kuk2)}2}, 

where  x'  =  /*  (d*x  +  n)  and  E  denotes  mathematical  expectation  [49].  Assume 
that  noise  n  is  a  random  field  with  zero  mean  and  is  independent  of  x,  the  expression 
of  £  can  be  evaluated  as: 

£   =   E[{f*d*x-x)2-2{f*d*x-x){f*  n)  +  (/  *  n)2] 
=  E  [((/  *  d  -  8)  *  xf]  +E[(f*  nf 

where  F,  D,  <3>xx  and  $„„  are  all  functions  of  (u^,^),  and  $xx  and  $„„  are  the  power 

density  spectrum  of  the  signal  and  noise  [49,  Section  2.10],  respectively.  Furthermore, 

the  integrand  can  be  factored  as: 

\FD  -  1|2  $xx  +  |F|2  $nn 
=  (F-  1/D)  (F*  -  1/D*)  \D\2  $xx  +  FF*<bnn 
=  (FF*  -  F/D*  -  F*/D  +  1/|D|2)  \D\2  $xx  +  FF*$nn 
=  FF*  (\D\2  $xx  +  $nn)  -  (F/D*  +  F*/D)  \D\2  $XI  +  $xx 


FF*  -  (F/D*  +  F*/D) 

_  (      \D\9XX  V 
\\D\2<t>xx+<t>nn) 


P  


£>|2$xx  +  $„„)  +  <I> 
(|D|2$xx  +  $ 

nn  -I 


F 


D'<S>XX 

\D\2<t>IX+<I>r 


(\D\2<S>    +$   1  lDl2*"     +  $ 
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Since  t  he  hist  term  is  large  than  or  equal  to  zero,  the  optimal  deblurring  filter 
(Wiener  filter)  is  thus 

F(ux,U2)  =  —   2     -,  (6.3) 

\D{u)^UJ2)\   9xx(Ui,U2)  +  «„n(Wi,W2) 

which  can  be  rewritten  in  terms  of  the  signal-to-noise  power  ratio 

D*{uuuj2) 


F(u}\,U}2)  —   9~ 

\D[u)X,U)2)\    +  $nn(uJl,UJ2)/$xx(uJuU2) 

For  white  Gaussian  noise,  $nn(^i, w2)  =  No-  The  power  spectrum  $xx{ui,u2) 
characterizes  correlation  property  of  the  signal  x.  As  the  correlation  usually  falls  off 
as  the  sample  distant  increases,  $xx(ui,  u2)  usually  is  a  lowpass  function.  If  the  signal 
has  no  correlation  among  its  adjacent  samples,  $xx(u)i,uj2)  =  S0,  and  the  signal-to- 
noise  power  ratio  equals  to  a  constant  S0/N0  =  1/R.  In  this  case,  the  Wiener  filter 
is  reduced  to 

%U2)  =  JM      .  (6.4) 

Comparing  with  the  gain-limited  inverse  filter  and  the  pseudo-inverse  filter, 
the  Wiener  filter  of  (6.4)  possesses  both  merits  of  continuity  and  noise  suppression. 

6.2.3    Wavelet- Vaguelette  Inversion 

The  combination  of  wavelet- vaguelette  decomposition  (WVD)  and  nonlinear 
shrinkage  of  WVD  coefficients  was  proposed  by  Donoho  [14]  for  linear  inversion  prob- 
lems involving  dilation-homogeneous  operators.  An  operator  T  with  such  property  is 
commutable  with  a  dilation  operator  Da  such  that 

TDa  =  aaDaT  , 

where  a  is  a  constant,  and  Da  is  the  dilation  operator  defined  as 

Da  [/(*)]  =  f{at)  . 

The  WVD  involves  three  sets  of  functions:  an  orthogonal  wavelet  basis 
{ipj,k{t)  =  2j/2*P(2jt  -  k)}  and  two  near-orthogonal  sets  {ujJe(t)  =  2j/2u{2jt  -  k)}  and 
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{vj,k(t)  =  2j/2v(2H  -  k)}.  The  three  sets  are  related  by  the  quasi-singular  value  re- 
lations 

r^j,fc  =  IjVjJe,  (6-5) 
TX*  =  7A*>  (6-6) 

where  quasi-singular  values  jj  depend  on  resolution  index  j  but  not  spatial  index  A;. 
Moreover,  the  two  sets  {%*(<)}  and  {vjtk(t)}  satisfy  biothogonality 

/oo 
Mj,*(*)Wm,n(')^  =  SJ.mh,n, 
-oo 

and  near-orthogonality 


2 

i,fc- 


Linking  everything  together,  the  reproducing  formula  of  the  WVD  is  expressed 

by 

=  S  trx'  7TViJfc(*)i 

where  [x,  y]  =  ff^,  x(i)y*(£)di  is  the  inner  product  of  x  and  y. 

Finally,  the  WVD  inversion  of  a  white  Gaussian  contaminated  measurement 
y(t)  =  Tx(t)  +  n(t)  may  be  achieved  via 

x'(t)  =  E^{b.M771}^(*)>  (6-7) 
j,* 

where  /C(j  is  the  nonlinear  shrinking  operator  defined  by  (4.4).  Notice  that  the 
threshold  tj  depends  only  on  resolution  j. 

6.2.4  Discussion 

The  Wiener  filter  provided  a  powerful  tool  for  linear  inversion  problems.  Given 
the  statistics  of  both  image  and  noise  processes,  one  can  design  a  optimal  filter. 
However,  Wiener  filter  belongs  to  the  linear  filtering  paradigm  and  works  in  the 
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frequency  domain.  The  linear  property  determines  that  it  has  to  trade-off  between 
smoothing  out  noise  and  sharpening  up  edges. 

In  the  opposite,  the  wavelet- vaguelette  approach  deployed  nonlinear  shrinkage 
to  combat  noise  in  the  time  domain.  It  is  characterized  by  three  components:  1) 
dyadic  analyzing  and  synthesizing  functions  determined  by  the  blurring  operator  T; 
2)  level-dependent  gain  factors  7"1;  and  3)  nonlinear  shrinkage  with  level-dependent 
thresholds  tj.  The  nonlinearity  enables  it  to  achieve  better  compromise  between 
denoising  and  sharpening.  However,  it  is  limited  to  dilation-homogeneous  operators 
only,  excluding  the  most  commonly  encountered  Gaussian  and  uniform  blurs  [14]. 
Moreover,  the  theory  was  developed  in  the  continuous-time  domain. 

The  limitation  of  the  WVD  inversion  lies  on  its  insistence  of  the  quasi-singular 
value  relations  (Eqs.  (6.5)  and  (6.6)).  For  a  convolution  operator  T  with  kernel  h(t), 
in  the  frequency  domain  (6.5)  is  equal  to 

H(u)^(2-jlo)  =  jjV{2-Jio).  (6.8) 
Let  u)'  =  2kuj  and  plug  it  into  (6.8),  we  have 

H(2ku)V{2-u-k)oj)  =  jjV(2-{j-k)u).  (6.9) 
In  the  other  hand,  directly  plugging  level  index  j  —  k  into  (6.8)  produces 

H{uj)y(2-{j-k)uj)  =  ^kV(2~^-k)u).  (6.10) 

By  comparing  (6.9)  and  (6.10),  we  conclude  that  it  must  be  true  that 

H(2ku)       7*  , 
TT.  ,    =  — —  (independent  of  ui)  , 
H(u)  7j_fc 

which  is  the  characteristic  of  dilation-homogeneous  operators  in  the  frequency  do- 
main. 
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6.3    Discrete-time  Overcomplete  Wavelet  Packet  Inversion 
6.3.1    One-dimensional  Formulation 

There  were  two  major  considerations  in  onr  searching  for  a  deblurring  algo- 
rithm. First,  it  should  be  available  in  the  discrete-time  domain.  Second,  it  should 
be  applicable  to  inhomogeneous  blurring  operators.  The  result  was  the  discrete-time 
overcomplete  wavelet  packet  representations  equipped  with  a  nonlinear  shrinking 
operator  and  a  gain  operator.  For  an  one-dimensional  observation 
y{k)  =  d(k)  *  x(k)  +  n(k),  our  deblurring  algorithm  can  be  expressed  by: 

00 

Wj{k)=   E  Vj{m)y{k-m),  (6.11) 

m=— oo 
A/-1  oo 

=  E  E  Kh  K(™)]  xJuj(k  -  m)>  (6-12) 

j'=0  m=  — oo 

where  {Aj|0<j<M-i}  is  the  gain  vector  and  K,tj  is  the  shrinking  operator.  Notice  that 
for  shrinking  operators  defined  by  (4.4)  it  is  easy  to  show 

JCtj  [wj(jn)}  Xj  =  Kfafo  [Wj(m)\j]  , 

which  means  that  the  order  of  gain  and  shrinking  operator  can  be  exchanged.  If  the 
shrinking  threshold  is  adaptive  to  magnitude  of  coefficients,  such  a  change  of  order 
would  not  alter  the  result.  We  prefer  the  order  of  "denoising"  and  then  "inversing"  as 
expressed  in  (6.12)  over  the  order  of  "inversing"  and  then  "denoising"  of  (6.7)  since 
we  believed  the  former  is  more  logical. 

The  idea  behind  our  method  is  the  same  as  in  the  WVD  inversion  (6.7).  That 
is  to  cut  off  noise  by  nonlinear  means  and  then  sharpen  the  survived  coefficients. 
Although  the  formula  (6.12)  looks  virtually  the  same  as  the  WVD  inversion,  it  has 
several  major  differences.  First,  it  uses  discrete-time  overcomplete  wavelet  packet 
representation,  which  guarantees  it  will  not  fall  on  the  anomaly  of  aliasing  enhance- 
ment. Second,  its  analyzing  and  synthesizing  functions  do  not  subject  to  the  dyadic 
bandwidth  constraint,  and  do  not  necessarily  satisfy  the  quasi-singular  value  relations 
(Eqs.  (6.5)  and  (6.6)). 


73 


To  complete  the  framework,  we  still  need  to  resolve  two  issues.  First,  for  a 
given  blur  with  frequency  response  D(u>),  how  can  we  determine  the  channel  config- 
uration (analyzing  and  synthesizing  filters)  as  well  as  the  gain  vector?  Second,  how 
should  we  choose  thresholds  tj  for  denoising? 

We  proposed  the  following  algorithm  for  the  channel  structure  and  the  gain 

vector: 

1)  Choose  a  gain-limited  inverse  filter  or  a  Wiener  filter  to  approximate. 

2)  Choose  the  prototype  filters  H(u),  H(u),  G(u)  and  G{u),  and  then  the  channel 
filters  of  the  Level  1  can  be  determined  by  Ci,0(^)  =  H(uj)H{uj)  and 
Ci,i{u)  =  G(u)G{u). 

3)  Call  Algorithm  4.2.2  (Minimum  Tree)  to  determine  both  the  tree  configuration 
and  gain  vector. 

However,  the  issue  of  determining  thresholds  tj  of  denoising  is  much  more 
difficult.  Thresholding  schemes  for  several  particular  applications  were  suggested  by 
Donoho  [15].  Generally  speaking,  knowledge  about  signals  and  noise  is  needed  to 
devise  a  thresholding  scheme.  By  treating  signal  and  noise  as  two  random  processes, 
the  correlation  property  may  be  used  to  characterize  them,  as  in  the  case  of  Wiener 
filtering.  For  wide-sense  stationary  white  Gaussian  noise  and  colored  Gaussian  signals, 
their  autocorrelation  functions  d>  and  power  spectrum  densities  $  may  be  written  as: 


Note  that  <3>XI(u;)  is  generally  different  from  the  power  density  spectrum  |X(u>)|2  of 
a  particular  realization  x  of  the  signal  process.  Under  these  assumptions,  a  possible 
thresholding  scheme  is: 


(6.13) 


(6.14) 
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where  Bj  is  the  bandwidth  of  the  channel  j,  and  uij  is  the  center  frequency  of  the 
channel. 

6.3.2    Two-dimensional  Extension 

The  one-dimensional  algorithm  of  the  overcomplete  wavelet  packet  inversion 
need  to  be  extended  to  two-dimension  in  order  to  deal  with  images  degraded  by 
symmetric  2D  blurs.  Formally,  the  2D  version  may  be  written  as 

A/-1  M-l 

x'{kuk2)  =J2  Y,  5m^c1,CJ[«;ci)c2(mi,m2)] ACl)C2uCl)CJ(fci  -m1,k2-m2). 

Cl=0  C2=0  mi  m2 

However,  an  efficient  2D  overcomplete  wavelet  packet  representation  is  much 
more  difficult  to  seek  than  in  the  ID  cases.  In  general,  the  2D  separable  overcomplete 
wavelet  packets  may  not  be  the  efficient  representation.  This  is  due  to  the  fact 
that  both  2D  modified  inverse  filters  and  Wiener  filters  are  non-separable  even  for 
separable  symmetric  2D  blurs.  For  example,  for  a  separable  Gaussian  blur: 

D(ux,uy)  =  Ae-^'+u2y^2, 

the  pseudo-inverse  filter  and  the  Wiener  filter  are 

l/AetA+tU*   ,  <  osjMA/R)  (615) 

0  ,  otherwise, 

and 

F(^u;y)  =  l/(l  +  R/Ae^+^°2). 

Both  filters  are  isotropic  and  can  be  efficiently  represented  by  donut-shaped  channel 
filters. 

In  practice,  a  reasonable  approximation  may  be  achieved  by  applying  the 
ID  algorithm  separately  to  rows  and  columns.  For  the  separable  Gaussian  blur,  a 
separable  pseudo-inverse  filter  may  be  constructed  as: 

F(ux,ujy)  =  Fi(u)x)Fi(ujy), 


F{UX,LOy) 
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where  the  ID  filter  Fi(lj)  is 
Fl{u) 


f  Ji/Ae"2/*2  ,|W|  <  a^HA/R) 
1  0  ,  otherwise. 


It  is  clear  that  the  rectangular  region  (\ujx\  <  a^/ln(.4/i?),  <  CT^ln(A/i?))  of  the 
above  2D  filter  includes  the  disk  of  yju^  +  to*  <  o\J\n(A/R)  of  the  (6.15).  By  choos- 
ing appropriate  parameter  R,  acceptable  results  may  also  be  achieved  by  approximat- 
ing 2D  Wiener  filters  with  ID  Wiener  filters.  For  uniform  blurs,  gain-limited  inverse 
filters  are  similar  to  Wiener  filters. 

The  major  problem  of  such  separable  approximations  is  they  may  boost  di- 
agonal frequency  components  significantly.  Figures  6.1  and  6.2  show  examples  of 
both  2D  non-separable  and  separable  filters.  The  emphasis  of  the  diagonal  frequency 
region  is  particularly  clear  in  the  case  of  Gaussian  blur. 

6.4    Experimental  Results 

To  quantitatively  compare  deblurring  results,  we  used  an  objective  measure, 
the  improvement  in  signal-to-noise  ratio  (ISNR),  defined  as  [3]: 

75^  =  20  logjlti},p=l,2, 

where  \\v\\p  is  Lp  norm  defined  as  [56]: 

/N-i  \  llv 


E  Hn)\P 


,n=0 


and  x,  y  and  x'  are  the  original,  degraded  and  deblurred  signals,  respectively. 

Moreover,  a  performance  upper  bound  is  useful  in  providing  a  gauge  for  com- 
parison. The  performance  of  an  oracle  Wiener  filter  may  be  considered  the  upper 
bound  for  Wiener  filters.  An  oracle  Wiener  filter  is  a  Wiener  filter  with  $nn(a>)  and 
$^(0;)  replaced  by  power  density  spectrums  |iV(a>)|2  and  lA'Xw)!2  of  the  particular 
realization  of  the  noise  and  the  original  signal.  The  filter  may  be  written  as: 

D*(u) 


|D(u;)|2  +  |aW/|A>) 


7G 


Row  1 


2 


3 


Figure  6.1:  Log  Frequency  Responses  of  Filters  for  a  Gaussian  Blur. 
Where: 

Row  1:  Blurring  filter 

Row  2:  2D  non-separable  filters,  gain-limited,  pseudo-inverse  and 
Wiener  filter  from  left  to  right 

Row  3:  2D  separable  filters,   gain-limited,   pseudo-inverse  and 
Wiener  filter  from  left  to  right 
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Figure  6.2:  Log  Frequency  Responses  of  Filters  for  an  Uniform  Blur. 


Row  1:  Blurring  filter 

Row  2:  2D  non-separable  filters,  gain-limited,  pseudo-inverse  and 
Wiener  filter  from  left  to  right 

Row  3:  2D  separable  filters,   gain-limited,   pseudo-inverse  and 
Wiener  filter  from  left  to  right 
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However,  it  is  not  necessarily  the  upper  bound  of  the  overcomplete  wavelet  packet 
inversion. 

6.4.1    One-dimensional  Signals 

One-dimensional  deblurring  examples  are  shown  in  Figure  6.3  and  6.4.  The 
ideal  signal  in  both  cases  is  a  perfect  multi-level  blocky  function,  which  is  very  chal- 
lenging and  revealing  for  any  deblurring  method.  For  overcomplete  wavelet  packet 
representations  used  in  both  examples,  Lemarie  filter  of  p=  1  was  used.  In  Figure  6.3, 
the  impulse  response  of  a  Gaussian  blur  was: 


Aexv[-{k/(aN/2))2},  0<k<N/2, 
Aexp[-((N  -  k)/{aN/2))2},   N/2<k  <N, 


with  a  =  0.03  and  ,4  =  1.0/  £fc  f(k).  The  impulse  response  of  the  uniform  blur  was 
/  =  [1, 1, 1, 1, 1, 1,  l]/7.  In  both  cases,  white  noise  with  a2n  =  0.01  was  added  such  that 
the  ratio  of  peak  signal  power  to  average  noise  power  was  SNR  =  101og(a;^ai/o-^) 
=  21.6d6.  To  make  fair  comparisons,  the  same  Gaussian  signal  model  of  (6.13) 
with  a  =  1.0  and  ax  =  0.007  was  used  for  both  Wiener  filtering  and  the  overcom- 
plete wavelet  packet  inversion.  The  thresholding  scheme  of  (6.14)  with  a  =  0.25  and 
a  =  0.18  was  used  for  the  Gaussian  and  uniform  blur,  respectively.  The  gain-limited 
inverse  filter  (6.1)  with  R  =  0.04  and  R  =  0.05  was  used  for  the  Gaussian  and  uni- 
form blur,  respectively.  For  the  minimum  tree  algorithm,  the  maximum  depth  was 
set  to  Level  6.  The  algorithm  found  a  configuration  of  twenty-five  channels  with 
e  =  5.7%  in  the  case  of  uniform  blur,  and  twelve  channels  with  e  =  0.004%  in  the 
case  of  Gaussian  blur. 

Table  6.1  compares  deblurring  performance  in  terms  of  I  SNR.  For  the  two 
particular  blurring  sources  and  the  signal,  it  shows  that  overall  performance  of  the 
overcomplete  wavelet  packet  inversion  is  about  the  same  as  Wiener  filtering  by  the 
measure  of  ISNR2.  and  it  is  better  than  Wiener  filtering  by  ISNRi. 
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Row  1 


Figure  6.3:  One-dimensional  Deblurrings  of  a  Gaussian  Blur. 

Where: 

Row  1:  The  original  signal 

Row  2:  The  blur  filter  and  the  degraded  signal 

Row  3:  The  oracle  Wiener  filter  and  the  deblurred  signal 

Row  4:  The  Wiener  filter  and  the  deblurred  signal 

Row  5:  The  gain-limited  inverse  filter  and  the  deblurred  signal 

Row  6:  The  overall  frequency  response  of  the  overcomplete  wavelet 
packet  inversion  and  the  deblurred  signal 

Row  7:  Channels  of  the  overcomplete  wavelet  packet  representation 

Note:  The  left  column  are  frequency  responses.  The  central  column 

are  signals  of  full-length.  The  right  column  are  zoom-in  of  the  feature 

region,  where  deblurring  results  were  overlaid  with  the  degraded  signal 

of  doted  line. 
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Figure  6.4:  One-dimensional  Deblurrings  of  an  Uniform  Blur. 

Where: 

Row  1:  The  original  signal 

Row  2:  The  blur  filter  and  the  degraded  signal 

Row  3:  The  oracle  Wiener  filter  and  the  deblurred  signal 

Row  4:  The  Wiener  filter  and  the  deblurred  signal 

Row  5:  The  gain-limited  inverse  filter  and  the  deblurred  signal 

Row  6:  The  overall  frequency  response  of  the  overcomplete  wavelet 
packet  inversion  and  the  deblurred  signal 

Row  7:  Channels  of  the  overcomplete  wavelet  packet  representation 

Note:  The  left  column  are  frequency  responses.  The  central  column 

are  signals  of  full-length.  The  right  column  are  zoom-in  of  the  feature 

region,  where  deblurring  results  were  overlaid  with  the  degraded  signal 

of  doted  line. 
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Table  6.1:  Performance  of  ID  Deblurring  Examples 


Blur  Type 

ISNR2  (DB) 

ISNRr  (DB) 

OWF 

\VF 

GLIF 

OWPI 

OWF 

WF 

GLIF 

OWPI 

Gaussian 

2.37 

1.88 

-9.04 

1.90 

-0.75 

-1.45 

-15.87 

1.18 

Uniform 

5.17 

3.90 

-4.32 

3.62 

-1.00 

-2.40 

-11.80 

0.00 

Note: 

OWF=Oracle  Wiener  Filter 
WF=Wiener  Filter 
GLIF=Gain-limited  Inverse  Filter 
OWPI=Overcomplete  Wavelet  Packet  Inversion 


6.4.2    Two-dimensional  Images 

Examples  of  image  deblurring  using  the  standard  Lena  image  were  included 
in  Figures  6.5  and  6.6.  Separable  blurring  sources  were  used.  The  same  ID  blurring 
filters  used  in  the  ID  examples  were  applied  row- wise  and  column- wise  on  the  original 
Lena  image.  A  statistical  signal  model  of  (6.13)  with  a  =  3.5  and  ox  =  0.008  was 
used  for  both  Wiener  filtering  and  the  overcomplete  wavelet  packet  inversion.  In  both 
cases,  white  Gaussian  noise  with  d2n  —  0.5  was  added  such  that  SNR  =  44.3DR  The 
thresholding  scheme  of  (6.14)  with  a  =  25.2  was  used  for  both  cases.  The  shrinking 
operator  was  applied  on  the  2D  components.  The  gain-limited  inverse  filter  (6.1)  with 
R  =  0.08  was  used  for  both  comparison  and  overcomplete  wavelet  packet  inversion.  A 
separable  overcomplete  wavelet  packet  inversion  algorithm  was  used  to  approximate 
a  2D  overcomplete  wavelet  packet  inversion.  For  the  case  of  Gaussian  blur,  the 
approximation  used  11x11  channels  with  e  =  0.005%  (ID).  For  the  case  of  uniform 
blur,  the  approximation  used  25  x  25  channels  with  e  =  3.9%  (ID).  The  deblurred 
images  were  all  linear  converted  to  the  same  intensity  range  of  the  original  image  for 
the  calculation  of  ISNR  indexes. 

Table  6.2  compares  deblurring  performance  in  terms  of  ISNR.  For  the  two 
particular  blurring  sources  and  the  Lena  image,  it  shows  that  overall  performance  of 
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Table  6.2:  Performance  of  2D  Deblurring  Examples 


Blur  Type 

ISNR2  (DB) 

ISNRi  (DB) 

OWF 

WF 

GLIF 

OWPI 

OWF 

WF 

GLIF 

OWPI 

Gaussian 

3.20 

3.17 

-14.04 

2.80 

2.96 

3.25 

-16.82 

2.69 

Uniform 

3.55 

1.78 

-14.78 

1.76 

L.46 

-0.60 

-18.05 

-0.64 

Note: 

OWF=Oracle  Wiener  Filter 
WF=Wiener  Filter 
GLIF=Gain-limited  Inverse  Filter 
OWPI=Overcomplete  Wavelet  Packet  Inversion 


the  overcomplete  wavelet  packet  inversion  is  about  the  same  as  Wiener  filtering  by 
both  measures  of  ISNR2  and  ISNRi. 

6.5    Summary  and  Discussion 

We  have  demonstrated  that  the  deblurring  problem  can  be  handled  using  an 
overcomplete  wavelet  packet  representation.  We  showed  that  inverse  filtering  in  the 
frequency  domain  and  denoising  in  the  time  domain  can  be  done  in  a  single  step.  In 
other  words,  for  this  application  the  gain  and  shrinking  operators  can  be  combined  as 
a  single  operation  on  the  spatial-frequency  representation.  Moreover,  our  approach 
may  be  seen  as  an  extension  of  the  wavelet- vaguelette  inversion  into  both  the  discrete- 
time  and  broader  problem  domains. 

Although  the  one  step  inversion  and  denoising  may  be  conceptually  appealing, 
its  limitation  is  also  apparent.  Since  nonlinear  shrinkage  works  only  within  areas 
where  the  magnitude  of  signal  components  is  stronger  than  that  of  noise  components, 
the  best  representation  for  denoising  is  to  maximize  local  signal-to-noise  ratio.  For 
example,  if  a  sharp  edge  is  decomposed  into  many  channels,  its  components  at  some 
channels  may  be  indistinguishable  from  noise.  Therefore,  the  best  representations 
for  approximating  the  inverse  filter  and  denoising  are  generally  different.  Having 
only  one  choice  of  representation,  the  approach  has  to  limit  its  application  to  certain 
signals  or  blurring  sources. 


Figure  6.5:  Deblurring  a  Gaussian-blurred  Lena  Image. 
Legend: 

(a)  The  original  256  x  256  image 

(b)  Degraded  image  with  SNR=44.3  db 

(c)  Deblurred  by  the  oracle  Wiener  filter 

(d)  Deblurred  by  a  Wiener  filter 

(e)  Deblurred  by  a  gain-limited  inverse  filter 

(f)  Deblurred  by  a  separable  overcomplete  wavelet  packet 
inversion 


(a)  (b) 


(c)  (d) 


(e)  (f) 

Figure  6.6:  Deblurring  an  Uniform-blurred  Lena  Image. 
Legend: 

(a)  The  original  256  x  256  image 

(b)  Degraded  image  with  SNR=44.3  db 

(c)  Deblurred  by  the  oracle  Wiener  filter 

(d)  Deblurred  by  a  Wiener  filter 

(e)  Deblurred  by  a  gain-limited  inverse  filter 

(f)  Deblurred  by  a  separable  overcomplete  wavelet  packet 
inversion 
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This  limitation  may  be  the  reason  that  our  preliminary  experimental  results 
showed  roughly  the  same  performance  as  the  much  simpler  approach  of  the  Wiener 
filtering.  For  the  Gaussian  and  uniform  blurring  sources,  some  channels  with  very 
narrow  bandwidth  are  needed  to  achieve  satisfactory  approximation  using  the  gain 
vector.  For  those  channels,  shrinking  operation  is  very  close  to  multiplication. 


CHAPTER  7 
CONCLUSION 

We  examined  both  classes  of  Fourier-based  and  wavelet-based  representations. 
We  showed  that  orthogonal  wavelet  transforms  lack  translation  invariance  and  the 
efficacy  of  representing  convolution  operators.  In  particular,  we  demonstrated  that 
aliasing  distortion  due  to  critical  subsampling  limits  its  applications  in  image  en- 
hancement and  restoration.  On  the  other  hand,  overcomplete  wavelet  representations, 
although  not  as  efficient  as  non-redundant  representations,  avoid  those  shortcomings. 
Moreover,  overcomplete  wavelet  representations  can  be  conveniently  analyzed  using 
linear  time  invariant  system  theory.  In  this  thesis,  we  evaluated  their  capability 
of  time-frequency  localization  using  the  uncertainty  factor.  We  also  identified  gain, 
shrinking  and  envelope  operators  to  approximate  convolution  operators,  suppress 
noise  and  extract  power  density  distributions.  We  investigated  the  criteria  of  mini- 
mum mean  square  error  and  designed  an  optimization  algorithm.  We  extended  the 
wavelet  shrinkage  to  overcomplete  representations.  We  also  presented  two  envelope 
detection  algorithms  and  compared  their  performances. 

The  framework  of  the  overcomplete  wavelet  representation  was  applied  to 
two  difficult  problems:  segmentation  of  textured  images  and  image  deblurring.  We 
demonstrated  that  channel  envelopes  as  feature  vector  produced  satisfactory  results 
on  both  natural  and  synthetic  textures.  We  showed  the  gain  and  shrinking  operators 
may  be  used  for  image  deblurring  and  pointed  out  its  limitations. 

The  issues  of  time-frequency  representations,  pattern  recognition  and  noise 
suppression  are  very  fundamental  ones.  Our  knowledge  about  them  is  going  to  be  ac- 
cumulative and  evolving.  This  thesis  increases  our  level  of  understanding  and  points 
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to  many  opportunities  for  future  research.  The  issue  of  time-frequency  localization 
was  touched  upon  but  not  fully  explored.  Although  the  uncertainty  factor  was  used 
for  evaluation  of  individual  channels,  we  lack  a  criteria  of  optimality  on  a  representa- 
tion as  a  whole.  Once  such  a  criteria  is  found,  we  may  able  to  search  for  the  optimal 
representation  of  textured  images,  which  was  not  explored  here.  Denoising  is  cer- 
tainly another  front.  The  optimal  representation  and  adaptive  selection  of  thresholds 
are  going  to  be  intertwined  and  challenging. 


APPENDIX  A 
PROOFS  RELATED  TO  DISCRETE-TIME  WAVELETS 

Some  basic  proofs  related  to  (bi)orthogonal  discrete-time  wavelets  were  in- 
cluded here  in  order  to  make  this  thesis  self-sustained. 

Proposition  A. 0.1  The  necessary  condition  for  a  pair  of  discrete  sequences  to  be 
orthogonal  (as  defined  by  (A.l))  is  (A. 2): 

oo 

]T  a(2n  -  m)b(m  -  21)  =  \5{n  -  I)  (A.l) 

m=— oo 

A{e^)B{ejuJ)  +  A(ej{uJ+n))B(e^+7!))  =  2A.  (A.2) 

where  A  is  a  constant,  and  A(ejoJ)  and  B{ejuJ)  are  discrete-time  Fourier  transform  of 
a(n)  andb(n),  respectively. 

Proof:  By  changing  variable  k  =  m-  21,  the  left  side  of  (A.l)  is  changed  to 
T,T=-oo  a  [2(n  -  l)  ~  k]  Hk)-  Then  it  is  equivalent  to  prove  that 

£  a(2l  -  k)b(k)  =  \S(l). 

k—-oo 

Notice  that  the  left  side  of  the  above  is  a  2-fold  decimation  of  the  series 

oo 

p(l)  =  £  a(l-  k)b(k). 

k=—oo 

Since  p(l)  is  the  convolution  of  a(l)  and  b(l),  its  discrete-time  Fourier  transform 
is  known  as 

P{e>u)  =  A{eju})B(e^). 
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Using  the  result  on  M-fold  decimator  proved  in  [63].  Fourier  transform  of  p(2l) 
is  equal  to 

l-  [P(e^2)  +  P{e«M)'2)\  , 
and  therefore  the  Fourier  transform  of  (A.l)  is 

P(e^/2)  +  p(eJ(^)/2)  =  2A. 

Finally,  by  substituting  P(eJUJ)  =  A{e?w)B{e?w)  back,  we  proved 
A{eJ"/2)B(ejuJ?2)  +  A(^w+2ir)/2)B{e>^+2w)/2)  =  2A  ■ 
Proposition  A.0.2  The  necessary  condition  for  (A. 3)  is  (A. 4): 

53  a(2/  -  m)b{n  -21)+  £  c(2/  -  m)d(n  -  2/)  =  8{m  -  n)  (A.3) 


/=— oo 


!=— oo 


A(eJ'w)B(ejw)  +C(eJ'w)D(eJ'w)  =  2. 
A(eJ'w)B(eJ'(w+,r))  +  C'(eJ'w)D(eJ'('J,+,r))  =  0. 


(A.4) 


Proof:  Consider  both  sides  of  (A.3)  to  be  two-dimensional  series  of  index  (m,n), 
and  apply  two-dimensional  Fourier  transform  to  the  both  sides. 
Let  p(m,  n)  =  E^-oo  a(21  ~  m)b{n  -  21),  then 


00  oo 


P(eJ'w,,e,'w»)   =     53     E  p(m,n)e-jmuJxe 


m=— oo  n=— oo 


=  Y. 


l=-oo  Lm=-oo 

00 

e 

/=— oo 
oo 


53  6(n  -  2/)e-jnw» 


5]  a(2J  -  m)e-,TBW" 

00 

^  a(ra)ejmu;* 

l=  — 00 

=     53  2tt5(2o;x  +  2a;„  +  2fc7r)i4(e-iw*)B(ejh'»') 


^  0-j7(2wx+2w1/) 


A;=  — oo 
oo 


t=  — oc 
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Therefore,  the  2D  Fourier  transform  of  the  left  side  is 

£  7T<y(u;x  +  a;y  +  for)  [.4(e-J^)^(^)  +  C(e-^)D(e^)j 


fc=  — 00 

and  that  of  the  right  side  is 

oo         oo  oo 


Y     Y  ${m-n)e-jmj*e-jnuJ"  =  Y  e-M"*+*>») 

m=  —  oo  n=—oo 


Y  2it6(ujx  +  u!y  +  2kir). 

k=-oo 


Thus,  2D  Fourier  transform  of  (A. 3)  is 

Y  5(ux  +  tuy  +  Attt)  [A{e-jw*)B(e>uy)  +  C{e-jw')D{e>^)\ 

fc=— oo 

oo 

=   Y  25(iox+u)y  +  2k7r). 

k=  —  oc 

By  combining  terms  with  even  and  odd  index  of  k,  the  above  equation  can  be 
rewritten  as 

oo 

Y  6(ux  +uy  +  2/ctt)  [X{e?u',e?u»)  -  2 

k=—oo 
oo 

+       Y.  5(ux  +  uy  +  (2k  +  l)ir)X(e?w*,e>u»)  =  0. 

k=  —  oo 

where 

X(e?w*,f?u*)  =  A{e-iu*)B{e?w*)  +  C(e-jW)D(eJ'^). 

Using  the  properties  of  the  5  function,  the  coefficients  of  the  above  8  series 
must  satisfy 

X(e-i{f*7k*\etu)  =  X(e-jw,ejuJ)  =  2, 
X(e-j{uJ+{2h+l)n\ejul)  =  X{e-^+7I\ejuJ)  =  0. 

Finally,  by  substituting  Xie^1^  e^Wy)  back,  we  proved 

A(ejuJ)B{ejuJ)  +  C(ejuJ)D(ejuJ)   =  2, 
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f2)-|C(^)[-^2)-|PM 


A(w)B(2w)  -(|4 


f4)-C(2u;)£>M 


(a) 


(b) 


Figure  A.l:  Equivalent  Structures  for  (a)  Convolution  and  Decimation  and  (b)  Ex- 
pansion and  Convolution. 

Proposition  A.0.3  Two  cascaded  stages  of  convolution  and  decimation  (expansion) 
are  equivalent  to  a  single  stage  convolution  and  decimation  ( expansion)  as  illustrated 
in  Figure  A.l. 

Proof:  Assume  input  x(n)  and  output  y(n),  and  a(n)  and  b(n)  are  the  inverse  DFT 
of  A(u>)  and  B(u). 
For  (a), 


y(n)    =  E 


00 

m=  — oo 


Y,f=_0Ox{k)a{2m  -  k)\  b(2n  -  m) 
=   ££-00  *(*)  [E^=-oo  fl(2m  -  k)b(2n  -  m)] 
=   ESUo  *(*)  [ES-oc  «(4n  -  k  -  206(0] 
=  E£_ooX(fc)p(4n-*) 

where  p(n)  =  a(n  -  21)6(0.  or  ^V")  =  A(e^)B(e^). 

Similarly,  for  (b), 


y(n) 


=  E^=-oo  [E2L-00  *(*)c(m  -  2*)j  d(n  -  2m) 

=  Er=-oc  *(*)  [E^=-oo  c(m  -  2k)d(n  -  2m) 

=  YSL-oo  *(*)  [E/^-oo  c(0d(n  -  4*  -  2/) 

=  E£_oo*(*)9(n-4*) 


where  g(n)  =  Ei=_occ(Orf(™  -  2/),  or  Q(e*")  =  C{ej2uJ)D(ej") 
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vk{n) 


Ufc(n) 


vc(n) 


g(n)  -(J2)   i'c,i(n) 


«c,l(n)  9(n)  -i 


L  A(n)  -Q2)  »c,2(n) 


«c,2(n)  (f2V-  h(n)  -1 


uc(n) 


Figure  A. 2:  Illustration  of  Adding  One  More  Channel  by  Splitting  a  Channel  into 
Two. 

Proposition  A. 0.4  The  basis  functions  for  a  binary  tree  structured  filter  bank  can 
be  expressed  as: 

k  =  0 


fc-i 


V0(e>")  =  G(e>w), 
fc-i 

Vk(e?u)  =  J]  H(e>*u)G(e>2k"),  Uk{e?")  =  J]  H(e^2'w)G(e'2kw),    I  <  k  <  M  —  2 

1=0  1=0 

A/-1  M-l 

Vk-i(^)  =  I]  iTfe*'2'1"),         £/A/-i(e>'w)  =  I]  H(ej2'u'l         k  =  M-  1. 
/=o  /=o 

u^ere  M  >  2. 

For  orthogonal  wavelet  basis  satisfying  condition  (2.11),  we  have 
Uk(en  =  Vk*(en,   or,  uk(n)  =  v*k(-n). 


Proof:  By  direct  application  of  the  Proposition  A. 0.3.  ■ 

Proposition  A. 0.5  Basis  functions  generated  by  a  general  binary  tree  with  filters 
hi(n)  and  gi(n)  satisfying  (2.9)  satisfy  the  following  biorthogonality: 


Em=-oo  uc(™  -  ncl)vk(nkn  -  m)  =6(c-  k)S(l  -  n) 


(A.5) 


Proof:  By  induction.  For  filter  bank  of  2  channels,  u0  —  gain),  ux  =  h0(n),  v0  =  g0{n) 
and  vi  =  h0(n).  In  this  case,  the  orthonomality  of  the  basis  functions  is  the  given 
property.  Now,  assume  (A.5)  is  true  for  up  to  M  channels.  To  construct  &M+1  chan- 
nel filter  bank,  we  split  an  arbitrary  channel  c  into  two  as  illustrated  in  Figure  A. 2. 
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Since  transform  coefficient  wc(n)  before  splitting  is 

we(n)  =  Zm=-ocx(m)vc(ncn  -  m), 
the  new  coefficient  wci(n)  and  wC2(n)  are 

Wd(n)   =  YZ-oo^SWn-l) 

=  E^-oc  *(m)  [E^-oc  Vc(nJ  ~  m)gc(2n  -  I) 
=   Y^=-rx,x(m)vcl{2ncn  -  m), 


where  vcl(n)  =  E/=_oo  VM  ~  ncl)gc(l). 
Similarly, 


vain)  =  E/=_oo  vc(n  -  ncl)hc(l), 

«cl(n)  =  E;=_oo  uc(n  ~  nJ)9c(l), 

u<a(n)  =  E™-oouc(n  ~  ncl)hc(l). 

For  k  /  cl  and  A-  ^  c2 

Sm=-oo  uci(m  -  ncli)vk(nkn  -  m) 
=  Em=-oo  [£~-oo  «c(m  -  2nci  -  nel)gc(l)j  vk(nkn  -  m) 
=  EZ-oc  9c(l)  [Em=-oo  ucim  -  nc{2i  -  l))vk{nkn  -  m)] 
=  E(=_oo  9c{l)5{c  -  k)6(2i  -  I  -  n)  =  0. 

and  similarly, 


Em=-oo  Mc2(m  -  ncli)vk(nkn  -  m)  =  0. 

For  /c  =  cl, 

oo 

53  itci(m  -  ncii)vci(ncin  -  m) 

m=-oo 
oo 

=  E 


5^  uc(m  -  2nci  -  ncli)gc(lx) 


l\  ——oo 


OO  00 


E     E  9c{h)gc{h) 


h=— oo  i2=— oo 
oo  oo 


53  vc(2ncn  -  m  -  ncl2)9c(k) 

/2  =  -00 


53  «c(m  -  77c(2i  +  Zi))uc(nc(2n  -  Z2)  -  m) 


=   E     E  Pc(Zi)5c(W2i  +  /1-2n  +  /2) 


/i  =-oo  h=— oo 
oo 


=   E  5c(ii)Pc(2(n-z)-/i)^5(n-i). 

il=  — 00 


94 


Similarly,  we  can  prove 


£m=-oo  Ucl(™ 


ncii)vC2(ncin  —  m)  =  8{n  —  i), 
ncli)vc2{ncXn  -  m)  =  0, 
ncii)vci(ncin  —  m)  —  0  ■ 


We  have  thus  proved  that  basis  functions  of  the  new  filter  bank  constructed 
this  way  are  also  orthonomal.  Notice  that  orthonomality  hold  under  quite  general 
conditions: 

1)  Filters  gc,  hc,  gc  and  hc  may  be  different  with  those  used  in  other  channels,  as 
long  as  they  satisfy  condition  (2.9). 

2)  Splitting  can  be  done  arbitrarily. 


APPENDIX  B 

NUMERICAL  COMPUTATION  OF  BATTLE-LEMARIE  WAVELET 

Battle-Lemarie  wavelets  are  built  from  polynomial  spline  of  order  2p+  1  [42]. 
Although  they  are  not  compactly  supported,  they  do  possess  two  important  prop- 
erties: symmetric  and  orthogonal.  The  symmetric  property  is  especially  desirable 
for  some  image  processing  applications  such  as  edge  detection,  character  recognition 
[39],  and  texture  segmentation  [62],  but  is  mutual  exclusive  with  the  compact  support 
property  [13]. 

In  this  appendix,  we  presented  an  algorithm  to  compute  the  scaling  function, 
wavelet,  and  associated  quadrature  mirror  filters  of  Battle-Lemarie  wavelets. 

B.l    Functions  in  Finite  Terms 

The  Fourier  transform  of  the  scaling  function  (f>,  wavelet  xp,  and  the  associated 
quadrature  mirror  filter  h  may  be  written  as  [42]: 

=   /  (B.l) 


XT/  I  <\           e              S4p+4(2  +tt)  m  9x 

*P  W      =     ~2^2\   V  1    YV  (BJ> 

uj2p+2  \  £4p+4(u;)X;4p+4(2) 


Hp(u)  = 

where  p  is  a  positive  integer  and 


£4p+4(q;) 
24p+4£4p+4(2u;) 


(B.3) 


+  00  ^ 

^n{u)  =    Yl    7  1   oi,  V 
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Therefore,  the  computation  of  the  function  En(u)  is  the  key.  Fortunately,  we  have  a 
closed  form  for  n  —  2  [42], 


1  1_ 

,(u;  +  2A;7r)2~4sin2(!)' 


=  m 


and  a  closed  form  expression  of  the  function  £n(u;)  based  on  £2(u;): 

(_l)n  An-2 

Thus,  all  we  needed  is  an  algorithm  to  compute  derivatives  of  the  function  T,2(u), 
which  we  shall  derive  in  the  following: 

Proposition  B.l.l  For  any  integer  n>\,  the  nth  order  derivative  of  the  function 
E2(w)  has  the  form  of 

dn  v  f  .      gn(cosf)  m  , 

=  SF*(f)'  (B-4) 

where  #„(cos  |)  is  a  polynomial  of  order  n  and  satisfied  the  following  recursive  rela- 
tion: 

<7„+i(s)  =  ~2  -  x  )  -  -^9n(x)  *  x  (B.5) 

Proof:  By  mathematical  induction. 

For  n  =  1,  £S2(u;)  =  -0.25|$g^.  Thus,  0i(cosw/2)  =  -0.25 cos(u;/2). 

Assume  for  any  n  >  1,  (B.4)  is  true,  and  #n(cos  |)  is  a  polynomial  of  order  n.  Then, 

for  n  +  1 , 

d"+1E2(a;)  =  _o[_  /gn(cosf)\ 
du^1     ~  do;  Vsinn+2(f )  J 

£  (gn(cos  f ))  *  sin"+2  -g„(cos  % )  *  (n  +  2)  *  sin"+1(f )  *  cos(|)  *  i 


;„2n+4 


sin 


(I) 


ff  (^(cos  I))  *  sin  |-^*gn(cosf).cos| 


dcos  . 

sin<n+1>+2(f) 


-|^^(cos  I)*!1- cos2  f)  -  ^  *  ^(cos  I)  * cos  f 

sin(«+l)+2(|) 

gn+1(cos|) 

sin(n+D+2(|)' 
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where  gn+i(x)  =  —  \  d  9"^(1  -  x2)  -  ^^-gn{x)  *  x,  and  is  clearly  of  order  (n  4-  1) 


Proposition  B.1.2  The  coefficients  {an<k}  of  the  polynomial  gn(x)  satisfies  the  fol- 
lowing recursive  relation: 


ao,o  =  0.25 

Oo,i  =  0 

an,Q  —  —Q-§a>n-\,\ 

an,k  =  0.5  [(fc  -  n  -  2)an_1>*_i  -  (fc  +  ljan-i.fc+i] ,  l<k<n-2 

On,n-l  =  — l-5an-l,n-2 


a. 


n,n 


'O'n—l.n—l 


(B-6) 


Proof:  For  n  =  0,  a0,o  =  0.25,  a0,i  =  0.  Based  on  Proposition  B.l.l,  we  can  write 


n-l 


9n-l(x)  =  X]  On-l,^* 


A:=0 


and. 

Using  (B.5),  we  can  derive 
9n(x)    =  - 


9n{x)  =  X!  an,k^ 


(B.7) 


1  rf5„_i(x)         2  _  (n+  1) 

V1        ^    /  o 


dx 

n-l 


1  A:=0  Z  fc=0 


Y,  an-hkkxk  1  +  -      an-i,kkxk+1  - 


an-i,kX 


k+l 


"  A:=0 
^  n-2 


k=0 


k=0 


n-l 


an-i,M  *  (l  +  1)  *  xl  +  -  Y,{k  -  n  -  1)  *  a„_i,fcx 


k+l 


1=0 


k=0 


^  7i—2  n 

=   --a„_u  -  -  £  an_i,i+i  *  (/  +       +  -  £(/  -n-2)  *  an-U-iXl 
1  1  i=\  1  i=i 

1  1  n~2 

=    --an_M  +  -  Y,  [(I  ~  n  ~  2)  *  On-l.Z-1  l)On-l,/+l]  X 

*  1  1=1 

—  nan-l,n-2-C        ~  U>n-l,n-lX 


(B-8) 


Comparing  (B.8)  with  (B.7),  we  obtain  (B.6) 


Proposition  B.1.3  For  even  order  n  of  the  polynomial  gn(x),  only  terms  with  an 
even  index  are  nonzero.  For  odd  order  n  of  the  polynomial  gn(x),  only  odd  index 
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terms  are  nonzero.  (This  property  halves  the  number  of  calculations  needed  for  the 
coefficients.) 

Proof:  By  induction. 

For  n  =  1,  ali0  =  — 0.5ao,i  =  0,  a^i  =  —  a^o  =  —0.25. 

For  n  =  2,  a2,0  =  -0.5au  =  0.125,  a2,i  =  -1.5ai.0  =  0,  a2,2  =  -ai,i  =  0.25. 
Assume  for  n  -  1,  the  assertion  is  valid. 
For  even  n,  let  n  =  2m, 


02m,277i       —    ~~  o2m_ii2m_i 
For  2m  is  even,  2m  -  1  is  odd.  And  for  odd  /,  /  —  1  is  even.  Use  the  induction 

hypothesis,  all  the  odd  index  terms  are  zero. 

The  similar  proof  applied  to  odd  order  (n)  polynomial.  ■ 

Based  on  the  above  properties,  we  can  now  write  T,n(u)  as 


G2m,0  —    —0.5  0,2m-\,\ 

02m,i        =   0.5  [(I  -  2m  -  2)a2m_i)i-i 

02m,2m+l    =    —1-5  a2m_i)2m_2 


(/  +  l)a2m-i,;+i] 


£n(u>)  = 


(-I)"  gn-2(cOSf) 

(n-1)!     sin"  | 


(B.9) 


By  combining  (B.9)  and  (B.3),  we  obtain 


and 


g4p+2(sin  |) 
\|  94P+2{cosuj)' 


From  Eqs.  (B.9),  (B.l)  and  (B.2),  we  obtain, 


'sin(u;/2)r+'         (4p  +  3)! 


(u/2)  J  ^24P+V+2(cosf)' 


and  therefore, 


*p(0) 


(4p  +  3)! 


%(0)  =  0. 


\|  2^^4p+2(l) 
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Having  derived  analytical  expressions,  we  shall  be  able  to  calculate  these  three 
functions  of  any  order  at  any  frequency  point. 

B.2    Numerical  Computation  Using  FFT 

The  three  functions  that  we  investigated  in  the  last  section  are  continuous 
functions  in  the  frequency  domain.  To  compute  their  waveforms  in  the  time  domain, 
we  used  the  FFT  (Fast  Fourier  Transform)  approach. 

Hp(uj)  is  a  27r-periodic  function.  It  is  Discrete  Fourier  Transform  of  quadrature 
mirror  filter  h(n)  which  is  a  discrete  function  in  time  domain.  Recall  that  sampling 
in  unit  circle  of  frequency  domain  results  in  periodic  extension  of  h(n)  in  time  do- 
main. However,  if  the  sampling  period  is  small  enough,  aliasing  in  time  domain  is 
negligible  and  h(n)  obtained  will  be  accurate.  The  sampling  values  of  the  Hp(u)  can 
be  obtained: 


2p+2 

I  \ 

 0<k<N-l. 

g4p+2{cosfk) 

Functions  <3>(w)  and  are  Fourier  transforms  of  scaling  function  <j>(t)  and 
wavelet  ip(t),  respectively.  Notice  that  functions  <f>(t)  and  t/)(t)  are  continuous  func- 
tions in  time  domain.  By  numerical  method,  we  can  only  obtain  sampling  values 
in  time  domain.  Recall  that  DFT  (Discrete  Fourier  Transform)  of  sampling  values 
(j)(nT)  and  ip(nT)  are  related  to  the  Fourier  transform  of  continuous  functions  <j)(t) 
and        by  [49]: 

*D(e'»)  =  -  Y,  *(r--7jr),  'MO  =  -  £  *(=--=-), 

1  k=-oc      11  1  k=-oo       1  1 

where  T  is  the  sampling  period. 

If  the  sampling  frequency  is  high  enough  such  that  fs  —  f^fh,  where  //,  is 
the  highest  frequency  component  in  (j)(t)  or  ip(t),  aliasing  can  be  neglected.  In  this 
case,  we  can  write: 
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Table  B.l:  Coefficients  of  g24 (x)  (a)  and  Truncated 

(a)   


Impulse  Response  of  hi(n)  (b). 
(b) 


k 

a2A,k 

n 

hi(n) 

0 

3679416778537.75 

0 

0.54173576 

2 

475582801692177.0 

1 

0.30682964 

4 

8011347151626016.5 

2 

-0.03549798 

6 

40907559988277171.5625 

3 

-0.07780792 

8 

81471913154443867.5 

4 

0.02268462 

10 

69736261315895408.25 

5 

0.02974682 

12 

26157606919852821.0 

6 

-0.01214549 

14 

4123375115615239.5 

7 

-0.01271542 

16 

243456149817686.25 

8 

0.00614143 

18 

4245937488920.0 

9 

0.00579932 

20 

13213718716.5 

10 

-0.00307863 

22 

2097148.875 

11 

-0.00274529 

24 

0.25 

12 

0.00154624 

Note: 


hp(-n)  =  hp(n) 


To  use  FFT  to  compute  values  <p(nT)  and  ip(nT),  we  need  to  sample  functions 
$£)(e^)  and  *D(eJu;)  as: 


and, 


*D(N-k)  ,%<k<(N-l). 


2  ' 


1*  (2ik\ 

T      \  NT  j 
VD(N-k)    ,f  <  k<  (N-l). 


,0  <  k  <  f , 


Sampling  values  (j){nT)  and  ip(nT)  can  then  be  obtained  from  inverse  FFT  of 
the  <J>£)(A:)  and  respectively. 

On  the  issue  of  implementation,  we  notice  that  coefficients  of  the  polynomial 
gn(x)  is  growing  very  fast.  As  an  example,  the  coefficients  of  #24(2)  is  shown  in  Table 
B.l.  In  our  numerical  calculations,  64-bit  precision  was  used.  Notice  that  gn(x) 
appears  in  both  numerator  and  denominator  of  the  formulae  of  Hp(u).  This  suggests 
another  way  of  preventing  overflow  in  calculating  Hp(io)  by  using  k(n)gn(x),  where 
k(n)  is  an  appropriate  multiplier  and  should  be  embedded  into  the  iterative  process 
of  (B.6). 
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Figure  B.l:  Transition  Band  of  Filters  Hp{uj)  with  p  =  1,3,9. 


As  the  order  p  of  the  filter  increases,  the  filter  Hp(u>)  is  closer  to  the  ideal 
halfband  filter.  Figure  B.l  showed  transition  band  of  filters  with  p  -  1, 3, 9. 

Two  MATLAB  m-files  for  calculating  QMF  of  the  B-L  wavelets  are  listed  in 
the  following. 

function  [H,G]  =  Lfilter(p,N) 

%  This  function  build  a  Quadrature  Mirror  Filter  pair  H,  G 
%  H:  low-pass  filter;  G:  high-pass  filter 
%  p:  2p  +  1  is  the  order  of  Lemarie  wavelet 
%  N  :  length  of  the  filters 

L  =  4*p  +  2; 

C  =  gcoef(L);      %  length  of  C  is  L+l 
H  =  zeros(l,N); 
G  =  zeros(l,  N); 

H{1)  =  sqrt{2);     %  at  w=0,  H=v/2 
H{N/2  +  1)  =  0.0;    %  at  w  =  tt  ,H=0 
dw  =  2*  pi/N; 
for  %  =  1  :  N/2  -  1 

omega  =  i  *  dw; 

u  =  0; 

v  =  0: 

for  k  =  1  :  2  :  L  +  1 

v  —  v  +  C(A:)  *  cos(omega/2)A(k  —  l); 

u  =  u  +  C(k)  *  cos(omega)  A(k  —  1); 
end: 

+  1)  =  s<?r£(2)*cos(ome#a/2)A(2*p+  2)*sqrt(v/u); 

end; 

if (JV/2  +  2  :  N)  =  H(N/2  :  -1  :  2); 
omega  =  [0  :  iV  —  1]  *  dw; 

G  =  -exp{-j  *  omega).  *  [H(N/2  +  1  :  N),  H{\  :  N/2)]; 
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function  [c]  =  gcoef(n) 

c  =  zeros(l,  n+1); 
t  =  zeros(l,  n  +  1); 
c(l)  =  0.25; 
c(2)  =  0; 
if(  n  >  8  ) 
k  =  n/4; 

else 

A:  =  n; 
end; 

for  i=l:l:n 

for  j  =  0  :  1 :  i 
if  j  ==  0 

*(1)  =  -0.5*c(2); 
elseif(  j>0  &  (j  <=  i-2)  ) 

t(j  +  1)  =  0-5  *((j  -  i  -  2)*c(j)  -  {j  +  l)*c(j  +  2)): 
elseif  j  ==  (i  —  1) 

t(j  +  l)  =  -1.5*c(*-l); 
else 

t(j  +  1)  =  -c(»); 

end; 
end; 

if(  i  >  A;  )     %  apply  multiplier 

c(l  :  t  +  1)  =  t(l  :  t  +  l)/(t  -  fc); 
else 

c(l :  t  +  1)  =  f(l  :  i  +  1); 
end: 
end: 


APPENDIX  C 


NUMERICAL  COMPUTATION  OF  THE  UNCERTAINTY  FACTOR 

Eq.  (3.7)  defines  the  uncertainty  factor  of  a  lowpass  filter  in  the  continuous 
frequency  domain.  In  order  to  compute  uncertainty  factors  of  filters  in  a  filter  bank, 
two  issues  need  to  be  resolved.  First,  the  definition  for  a  lowpass  filter  has  to  be 
extended  to  bandpass  and  highpass  filters.  Second,  the  formulas  have  to  be  modified 
for  computation  in  discrete  frequency  domain. 

The  general  definition  of  the  variance  of  the  frequency  distribution,  includ- 
ing bandpass  filters,  adopted  by  Liu  and  Akansu  [40]  is: 


The  drawback  of  applying  this  definition  to  bandpass  filters  is  that  for  real 
filters  \F(eju>)\  =  \F(e~^)\  and  thus  u;  =  0  and  is  not  the  central  frequency  of  the 
passband.  Therefore,  this  is  not  a  fair  comparison  between  lowpass  and  bandpass 
filters.  Moreover,  the  variance  will  also  depend  on  the  location  of  the  passband 
due  to  the  factor  (u  —  u;)2  =  a;2,  and  thus  is  not  fair  even  among  bandpass  filters. 

To  overcome  this  problem,  we  notice  that  the  idea  behind  the  variance  aw 
is  the  energy  distribution  within  a  complete  passband.  For  bandpass  filters,  a  single 
passband  in  the  positive  frequency  side  should  be  used.  We  thus  adopted  the  following 
definition: 


where 


C  \F(ei«)\Muj 


(C.l) 
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Figure  C.l:  Channel  Bandwidths  of  Overcomplete  Wavelet  Representations. 

C  u\F{&)\*du> 


u,  = 


_  JuJi_ 


Q  \F{eJ»)\*du  ' 


(C.2) 


where 


Wl  =  < 


for  lowpass  filters, 
for  bandpass  filters, 
for  highpass  filters. 

Notice  that  this  definition  is  good  only  for  filters  whose  passband  is  completely  con- 


-7T, 
0. 
0. 


7T, 
27T, 


tained  in  the  range  [uji,ujh],  although  sidebands  are  allowed  within  this  range.  Also 
notice  that  the  bandwidth  of  the  lowpass  band  is  twice  the  width  of  bandpass  bands 
by  the  construction  of  (3.3)  as  illustrated  in  Figure  C.l. 

In  discrete  frequency  domain,  Eqs.  (C.l)  and  (C.2)  should  be  modified  as: 

a    =    Et"(fe-A;c)2|F(fc)|2  /27TX2 
Ekl  \F(k)\2       \N)  ' 

where  N  is  the  number  of  sampling  points  in  the  [0,  27r],  and 


k. 


'  -N/2, 
0, 
0. 


N/2,        for  lowpass  filters, 
N/2,        for  bandpass  filters, 
N,        for  highpass  filters. 


Notice  that  \F(k)\  is  N-periodic  such  that  \F(k)\  =  \F(k  +  N)\,  and  therefore 
\F(-k)\  =  \F(N-k)\. 

Similarly,  due  to  the  periodic  nature  of  the  sequence  f(n)  in  the  time  domain, 
the  variance  on  should  be  modified  as 

E^ol(«-n)2|/(n  +  d)|2 


mm 

0<d<N 


n  = 


gai/(n+*oia 

ZZ=Qln\f(n  +  d)\> 
Z^\f(n  +  d)\*  ' 
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where  f(n  +  d)  is  the  original  f(n)  circular  shifted  by  d. 

A  MATLAB  function  uf actor ()  for  computing  the  uncertainty  factor  is  in- 
cluded in  the  following.  Also  included  is  a  MATLAB  function  lwuf()  for  building 
analyzing  filters  of  an  overcomplete  wavelet  packet  using  a  prototype  Lemarie  filter 
and  computing  uncertainty  factors  of  them,  where  the  function  Lf  titer ()  is  given  in 
Appendix  B. 

function  [u]=ufactor(H) 

%  H  is  the  Discrete  Fourier  Transform  of  a  discrete  filter  h(n) 
N  =  length{H); 
W  =  H.*conj(H); 
vmax  =  max(W); 

if  W(l)  >  0.9b*vmax  %  low-pass  filter 

F  =  \W{N/2  -1:N),  Wi\  :  N/2)]; 
elseif  W(N/2)  >  0.9b*vmax  %  high-pass  filter 

F  =  W; 

else  %  band-pass  filter 

F  =  W{\  :  N/2); 
end; 

L  =  length{F); 

fc=[0:  L-l\*  F'/sum(F); 

sf  =  (([0  :  L-l)  -  fc).  A  2)  *  F' /sum(F)*(2*pi/N)  A2; 

%  transform  to  time  domain 
h  =  real(ifft(H)); 
st  =  lelO; 

for  d  =  2  :  N- 1       %  circular  shift  to  find  the  minimum 
x  =  [h{d:  N),h(l  :  d-1)]; 
\Y  =  x.  A  2; 

tc=[0  :  N-l}*  W'/sum(W); 

m  =  (([0  :  N-l]  -  tc).  A  2)  *  W'/sum{W); 

if(  m  <  st  ) 

st  =  m; 
end: 
end: 

u  =  st  *  sf; 

function  []  =  lwuf(shell,p,nlevel) 

%  shell  :  =1,  autocorrelation  shell;  =0,  regular 
%  p  :  2p  +  1  is  the  order  of  Lemarie  wavelet 
%  nlevel:  number  of  levels 
N  =  2  A  nlevel; 
SIZE  =  256; 

omega  =  2*pi*[0  :  SIZE-l]/ SIZE; 
[H.G]  =  Lf  titer  {p,  SIZE); 
H  =  H/sqrt{2); 
G  =  G/sqrt{2); 
if  shell  ==  1 

H  =  H.*  conj{H); 


G  =  G.  *  conj{G); 
end: 

W  =  ones(N,  SIZE); 
HH  =  ones(l,  SIZE); 
GG  =  ones (1,  SIZE); 
d=l; 

for  72  =  1  '.  hIcvgI 

HH  =  H(rem([0  :  {SIZE-l)]*d,  SIZE)  +  1); 
GG  =  G{rem([0  :  (SIZE-l)]*d,SIZE)  +  l); 
for  j  =  2An  :  -2  :  1 
Irf  =  ceil{j/2); 
if  rem(pt,  2)  =  0 

Wj,:)  =  W(pt,:).*GG; 
W(j-1,:)  =  W(pt,:).  * 

Gls6 

W(j,:)  =  W(pt,:).*HH; 
W(j-l,:)=W(pt,:).*GG; 
end: 
end; 

of  =  d*2; 
end: 

clg; 

mi  =  1.1  *  max(max(abs(W (:, :)))); 
for  n  =  1  :  AT 

f  =  abS(W(n,:)); 
if  n  <=  TV/2 

p\ot(om.ega,  f); 
else 

plot  (ome^a,  /,' :'); 
end; 

uf  =  u/a(ior(VT(n, :)); 

str  =  sprint /('%.2/',  «/); 

text(((n- 1)+  0.1)  *pi/iV,  0.95  *mrr,  sir); 

if  n==l 

axis([0  pi  0  mi]); 

hold; 
end; 
end; 

hold  off; 

of  =  sprintf('p%dl%ds%d.eps',p,  nlevel,  shell); 
eval(['print  -deps  '  of]); 


APPENDIX  D 

BOUNDARY  TREATMENTS  OF  FINITE-LENGTH  SEQUENCES 
D.l    Periodic  Extensions 

For  convolutions  involving  finite-length  discrete  sequences,  we  face  the  bound- 
ary problem.  One  solution  is  to  construct  an  infinite-length  signal  from  the  finite- 
length  sequence.  For  a  finite-length  sequence  s(n)  of  length  N,  two  constructions  are 
commonly  used: 

1)  N-Periodic  extension.   An  infinite-length  signal  s(n)  of  period  N  is  con- 
structed as: 

s(n)-{  S{U)  ,0<n<N-l  \  =  ~{n  +  kN) 

S[n)~\s(n  +  N)    ,-JV<n<-l  (     s[n  +  Kiv). 

This  method  is  equivalent  to  a  N-point  FFT  implementation  of  convolution. 

Since  the  periodic  extension  generally  introduce  discontinuities  at  boundaries, 

the  output  signal  usually  exhibits  large  artifacts  near  the  boundaries. 

2)  2N-Periodic  mirror  extension.  An  infinite-length  signal  s(n)  of  period  2iV 

is  constructed  as: 

-/  x     f  s(n)  ,  0  <  n  <  N-  1  1  AT, 

s{n)  =  {    )'  '     ~         .        \=s(n  +  2kN). 

I  s{—n  —  1)   ,  —  N  <  n  <  —  1 

Since  the  new  signal  s(n)  is  continuous  at  boundaries,  the  undesirable  artifacts 

may  be  significantly  reduced. 

D.2    Closure  of  Symmetry  Under  Convolution 

In  order  to  define  symmetric/antisymmetric  sequences,  we  first  define  an  ex- 
pansion operator  |  as: 

tm  (s(n))  = 


s(n/m)  ,  n  =  0,  ±m,  ±2m, 
0  ,  otherwise. 
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and  denote  tm  (s(n))  (n).  Equivalently,  we  may  write: 

00 

s^m(n)  =  E  s(k)5(n  —  km).  (D.l) 

k=  —  oo 

We  then  define  symmetric/antisymmetric  sequences  as: 

Definition  D.2.1  For  a  discrete  sequence  s(n),  if  there  exist  integers  c  and  m  such 
that  Sfm(c  —  n)  =  s^m(c  +  n)  (s^m(c  —  n)  =  —s^m(c  +  n)),  the  sequence  is  said  to  be 
symmetric  (antisymmetric)  at  c/m,  and  C  =  c/m  is  said  to  be  the  symmetry  center. 
For  symmetry  sequences,  we  may  define  a  symmetry  index  S:  S  =  1  for  symmetric 
sequences,  and  S=  -1  for  antisymmetric  sequences. 

Symmetry  sequences  possess  some  special  properties.  We  will  prove  that  sym- 
metry sequences  is  closed  under  convolution.  We  first  prove  a  proposition. 

Proposition  D.2.2  The  expansion  operation  is  distributive  over  convolution: 

r(s*f)(n)  =r  (s(n)hr  (/(*)). 

Proof:  According  to  (s  *  f)(n)  =  T.T=-oo  s{k)f(n  ~  k)  and  (D-1)*  we  have  the  left 
side: 


00  00 


r(s*f)(n)   =    E    E  s(k)f(l-k)S(n-lm) 

l=  —  oo  fc=— oo 
oo  oo 

=    E  s(k)  E  f(l-k)S{n-lm) 

fr  =  -oc  l=— oo 

00  oo 


E  «(*)  E  mS{n-(i  +  k)m) 

k=  —  oo  i=—oo 


and  the  right  side: 


00  00 


r(s(n))*r(f(n))   =     E    E  s(l)5(k  -  Im)  mS(n-k-im) 

fe=— oo /=— oo  i=— oo 

00  00  oo 

=     E  s(l)  E  /(O  E  *(*  -  lm)S(n  -  k  -  im) 

/=  — oo         i=  —  oo         fc=— oo 
oo  oo 

=    E  s(0  E  f{i)^(n  —  lm  —  im)  m 

l=  —  00  t=  —  00 

We  now  prove  the  closure  property  on  convolution. 
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Theorem  D.2.3  Convolution  of  two  symmetry  sequences  s(n)  and  f(n)  results  in 
a  symmetry  sequence  y(n)  with  symmetry  index  Sy  =  Sx  ■  Sj  and  symmetry  center 

Cy    CX    +  Cf  . 

Proof:  Based  on  the  given  condition,  we  have: 

x-^m  (cx  —  ?i)  —  Sxx^m  {cx  ~\~  nj 
fVn  (cf  -  n)  =  Sffrm  (cf  +  n) 

Applying  Proposition  D.2.2,  we  have: 

00 

Vtm  icx  +  c/  —  n)   =     ^2  X|m  (k) ffm  (cx  +  Cf  —  n  —  k) 

k=—oo 
00 

=    ]L  xtm  (c*  ~~  0/tm  (c/  -  n  +  0 
i=— 00 
00 

=    ]T  5xar-fm  (cx  +  l)Sffrm  (cf  +  n  -  I) 

l=—oo 

OO 

=   SxSf  ^2  xtm(k)f^m(cx  +  Cf  +  n  -  k) 

k=  —  oo 

=  SxSfy^m(cx  +  Cf  +  n)  ■ 

Since  translation  of  index  n  will  not  actually  alter  a  discrete  sequence,  without 
loss  of  generality,  symmetry  sequences  may  be  classified  according  to  their  symmetry 
index  S  and  symmetry  center  C  as  below: 

•  Type-I  symmetric.  s(n)  =  s(—n).  The  symmetry  center  C  =  0. 

•  Type-II  symmetric.  s(n)  =  s(—n  -  1).  The  symmetry  center  C  =  -1/2 
(c  =  —  \,m  =  2).  The  2N-periodic  mirror  extended  sequences  introduced  pre- 
viously belong  to  this  class. 

•  Type-Ill  symmetric.  s(n)  =  s(-n  +  1).  The  symmetry  center  C  =  1/2 
(c  =  l,m  =  2). 
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•  Type-I  antisymmetric.  s(n)  =  -s(-n).  The  symmetry  center  C  =  0.  No- 
tice that  s(0)  =0  is  the  constraint. 

•  Type-II  antisymmetric.  s(n)  =  -s(-n  -  1).  The  symmetry  center  C  =  -1/2 
(c=  -l,m  =  2). 

•  Type-Ill  antisymmetric.  s(n)  =  -s(-n+  1).  The  symmetry  center  C  =  1/2 
(c=  l,m=  2). 

The  closure  property  of  symmetric  sequences  may  reduce  the  computation 
burden  significantly.  For  symmetric  filters  and  mirror  extension,  we  only  need  to 
carry  out  convolution  on  at  most  N  +  1  sample  points  instead  of  2N. 
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